Cone-beam reconstruction using backprojection of locally filtered projections and X-ray CT apparatus

ABSTRACT

Embodiments of the present invention include a method and apparatus for accurate cone beam reconstruction with source positions on a curve (or set of curves). The inversion formulas employed by embodiments of the method of the present invention are based on first backprojecting a simple derivative in the projection space and then applying a Hilbert transform inversion in the image space.

CROSS-REFERENCE TO RELATED APPLICATIONS

This nonprovisional patent application claims benefit and priority under35 U.S.C. § 119(e) of the filing of U.S. Provisional Patent ApplicationSer. No. 60/694,652 filed on Jun. 28, 2005, titled “CONE-BEAMRECONSTRUCTION,” the contents of which are incorporated herein byreference for all purposes.

BACKGROUND OF THE INVENTION

1. Field of the Invention

This invention relates generally to reconstruction of the derisityfunction of a three-dimensional object from a set of cone-beamprojections, such as from an X-ray source. More particularly, theinvention relates to methods for cone beam reconstruction usingbackprojection of locally filtered projections and an X-ray computedtomography (CT) apparatus incorporating the method.

2. Description of Related Art

Technology for X-ray detection in cone-beam (CB) geometry is rapidlyimproving and offers more and more potential for the construction ofrobust computed tomography (CT) systems for fast high-resolution volumeimaging. However, to optimally build such systems, the problem of CBimage reconstruction needs to be fully understood.

The successive works of H. Tuy, “An Inversion Formula for Cone-BeamReconstruction,” SIAM J. Appl. Math., No. 43, pp. 546-52, 1983, B. D.Smith, “Image Reconstruction from Cone-Beam Projections: Necessary andSufficient Conditions and Reconstruction Methods,” IEEE Trans. Med.Imag., Vol. MI-4, pp. 1425, 1985 and P. Grangeat, “MathematicalFramework of Cone-Beam 3D Reconstruction via the First Derivative of theRadon Transform,” in Mathematical Methods in Tomography, G. T. Herman,A. K, Louis, and F. Natterer, Eds. Berlin, Germany: Springer-Verlag,1991, Vol. 1497, Lecture Notes in Mathematics, pp. 66-97, have shownthat exact reconstruction at a given location is possible if every planepassing through that location intersects the trajectory of the X-raysource. However, this well-known result, which is fundamental andclearly represented a breakthrough, is not sufficient for most practicalapplications because its derivative assumed complete CB projections.When only part of the imaged object is illuminated at a given sourceposition, the CB projection is said to be incomplete or truncated. Exactand accurate reconstruction from truncated projections is morecomplicated than from complete projections. An overview of the problemof reconstructing from truncated projections may be found in R.Clackdoyle, M. Defrise and F. Noo, “Early Results on General Vertex Setsand Truncated Projections in Cone-Beam Tomography,” in ComputationalRadiology and Imaging: Therapy and Diagnostics, C. Brgers and F.Natterer, Eds. Berlin, Germany: Springer-Verlag, 1999, Vol. 110, IMAVolumes in Mathematics and Its Applications, pp. 113-135. Under someconditions, the problem of reconstructing from truncated projections mayeven be impossible, see e.g., F. Natterer, The Mathematics of CT,Philadelphia, Pa.: SIAM, 2001.

A general theory to handle CB data truncation remains elusive. Solutionshave been found only for particular measurement geometries. Many ofthese solutions were obtained using a clever handling of data redundancyin a CB filtered-backprojection (FBP) reconstruction framework. Seee.g., H. Kudo and T. Saito, “An Extended Completeness Condition forExact Cone-Beam Reconstruction and Its Applications,” in Conf. Rec. 1994IEEE Nuclear Science Symp. Medical Imaging Conf., Norfolk Va., 1995, pp.1710-14; F. Noo, R. Clackdoyle, and M. Defrise, “Direct Reconstructionof Cone-Beam Data Acquired with a Vertex Path Containing a Circle,” IEEETrans, Image Process., Vol. 7, No. 6, pp. 854-67, Jun. 1998; R. H.Johnson, H. Hu, S. T. Haworth, P. S. Cho, C. A. Dawson and J. H.Linehan, “Feldkamp and Circle and Line Cone-Beam Reconstruction for 3DMicro-CT of Vascular Networks,” Phys. Med. Biol. Vol. 43, pp. 929-40,1998; H. Kudo and T. Saito, “Fast and Stable Cone-Beam FilteredBackprojection Method for Nonplanar Orbits,” Phys. Med. Biol., Vol. 43,pp. 747-60, 1998; H. Kudo, F. Noo and M. Defrise, “Quasi-Exact FilteredBackprojection Algorithm for Long-Object Problem in Helical Cone-BeamTomography,” IEEE Trans, Med. Imag., Vil. 19, No. 9, pp. 902-21, Sep.2000; G. Lauritsch, K. C. Tam, K. Sourbelle and S. Schaller, “ExactLocal Region-of-Interest Reconstruction in Spiral Cone-BeamFiltered-Backprojection CT: Numerical Implementation and First ImageResults,” in Proc. SPIE Medical Imaging Conf. (Image Processing), Vol.3979, 2000, pp. 520-32; K. C. Tam, G. Lauritsch and K. Sourbelle,“Filtering Point Spread Function in Backprojection Cone-Beam CT and ItsApplications in Long Object Imaging,” Phys, Med. Biol., Vol. 47, pp.2685-703, 2002; A. Katesevich, “A General Scheme for ConstructingInversion Algorithms for Cone-Beam CT,” I.J.M.M.S., Vol. 21, pp.1305-21, 2003; G. H. Chen, “An Alternative Derivation of Katsevich'sCone-Beam Reconstruction Formula,” Med. Phys., Vol. 30, No. 12, pp.3217-26, 2003; J. D. Pack, F. Noo and H. Kudo, “Investigation of SaddleTrajectories for Cardiac CT Imaging in Cone-Beam Geometry,” Phys. Med.Biol., Vol. 49, No. 11, pp. 2317-36; and H. Kudo, F. Noo and M. Defrise,“Cone-Beam Filtered-Backprojection Algorithm for Truncated HelicalData,” Phys. Med. Biol., Vol. 43, pp. 2885-909, 1998.

Other solutions for CB reconstructing from truncated projections forparticular measurement geometries were obtained using the formula posedby Grangeat or its truncated version, see H. Kudo, F. Foo, and M.Defrise, “Cone-Beam Filtered-Backprojection Algorithm for TruncatedHelical Data,” Phys. Med. Biol., Vol. 43, pp. 2885-909, 1998, incombination with properties of the three-dimensional (3-D) radontransform. See, e.g., H. Kudo and T. Saito, “Extended Cone-BeamReconstruction Using Radon Transform,” in Conf. Rec. 1996 IEEE NuclearScience Symp. Medical Imaging Conf., Anaheim, Calif., 1997, pp. 1693-97;K. C. Tam, S. Samarasekera and F. Sauer, “Exact Cone-Beam CT with aSpiral Scan,” Phys. Med. Biol, Vol. 43, pp. 1015-24, 1998; S. Schaller,F. Noo, F. Sauer, K. C. Tam, G. Lauritsch and T. Flohr, “Exact RadonRebinning Algorithm for the Long Object Problem in Helical Cone-BeamCT,”IEEE Trans, Med. Imag., Vol. 19, No. 5, pp. 361-75, May 2000; and X.Tang and R. Ning, “A Cone Beam Filtered Backprojection (CB-FBP)Reconstruction Algorithm for a Circle-Plus-Two-Arc Orbit,” Med. Phys.,Vol. 28, No. 6, pp. 1042-55, 2001.

Many advances in CB reconstruction have been made recently inreconstruction methods for use in helical CB tomography (HCBT). Two ofthe most interesting results achieved in the HCBT context are disclosedin A. I. Katsevich, “An Improved Exact Filtered Backprojection Algorithmfor Spiral Computed Tomography,” Adv. Appl. Math., Vol. 32, No. 4, pp.681-97, 2004 and Y. Zou and X. Pan, “Exact Image Reconstruction on PILines from Minimum Data in Helical Cone-Beam CT,” Phys. Med. Biol., Vol.49, pp. 941-59, 2004.

Katsevich showed that exact and accurate FBP reconstruction can beachieved with a simple one-dimensional (1-D) Hilbert transform of aderivative of the projection, using a minimum overscan and just slightlymore than the data in the Tam-Danielsson (TD) window, see e.g., K. C.Tam, S. Samarasekera and F. Sauer, “Exact Cone-Beam CT with a SpiralScan,” Phys. Med. Biol, Vol. 43, pp. 1015-24, 1998; and P. E.Danielsson, P. Edholm, J. Eriksson and M. Magnusson Seger, “Toward ExactReconstruction for Helical Cone-Beam Scanning of Long Objects: A NewDetector Arrangement and a New Completeness Condition,” in Proc. 1997Meeting Fully 3D Image Reconstruction in Radiology and Nuclear Medicine,D. W. Townsend and P. E. Kinahan, Eds., Pittsburgh, Pa., 1997, pp.141-44. Katsevich has also generalized his formula to general CBtomography and showed incidentally that his formula belongs to thefamily of methods that can be obtained using a clever handling of dataredundancy in a CB-FBP reconstruction framework, see e.g., A.Katesevich, “A General Scheme for Constructing Inversion Algorithms forCone-Beam CT,” I.J.M.M.S., Vol. 21, pp. 1305-21, 2003.

Zou and Pan investigated the effect of skipping the 1-D Hilberttransform in the reconstruction steps of Katsevich's formula. Theyargued that by so doing the outcome of the backprojection on any π-lineis the Hilbert transform of the values of the density function, ƒ, onthis π-line, and they devised from this argument and information on thesupport of ƒ an accurate algorithm that has the same overscan andefficiency as Katsevich's formula but uses only the data in the TDwindow. Recall that a π-line is any line segment that connects twopoints of the helix separated by less than one helix turn, see e.g., P.E. Danielsson, P. Edholm, J. Eriksson and M. Magnusson Seger, “TowardExact Reconstruction for Helical Cone-Beam Scanning of Long Objects. ANew Detector Arrangement and a New Completeness Condition,” in Proc.1997 Meeting Fully 3D Image Reconstruction in Radiology and NuclearMedicine, D. W. Townsend and P. E. Kinahan, Eds., Pittsburgh, Pa., 1997,pp. 141-44.

None of these conventional approaches appears to achieve CBreconstruction on various measured lines using virtually arbitrarysource trajectories. A measured line is any line that contains a sourceposition and is part of the measurements. Thus, it would be highlyadvantageous to provide a method for CB reconstruction usingbackprojection of locally filtered projections and an X-ray CT apparatusincorporating such a method.

SUMMARY OF THE INVENTION

Embodiments of the present invention achieve reconstruction on variousmeasured lines (M-lines) using (almost) arbitrary source trajectories.Embodiments of the present invention include applying the processingsteps of the well-known algorithm disclosed in L. A. Feldkamp, L. C.Davis and J. W. Kress, “Practical Cone-Beam Algorithm,” J. Opt. Soc. Am.A, Vol. 1, pp. 612-19, 1984, to segments of source trajectories using asimple derivative along the detector rows instead of the ramp filter andobtain a portion of the Hilbert transform of ƒ on various measuredlines. Embodiments of the present invention further include achievingthe reconstruction of ƒ on these measured lines by obtaining a 1-Dfunction from its Hilbert transform when the latter is only known overthe region where the 1-D function is nonzero.

An embodiment of a method to obtain ƒ on line segments that connect twosource positions is disclosed according to the present invention. Suchline segments are referred to herein as R-lines, which is a short-handnotation for redundantly measured lines. When the X-ray sourcetrajectory is a helix and the extremities of the R-line are separated byless than one helix turn, the R-line is a π-line. In general, there isno guarantee for a voxel to belong to an R-line, so the methods of thepresent invention do not allow reconstruction at arbitrary locations.The locations where reconstruction is achievable is dependent on thesource trajectory and the extent of the detector. Fortunately, a largenumber of data acquisition geometries have their field-of-view (FOV)entirely covered by R-lines, see, e.g., Danielsson et al. for the helixand Pack et al. for the class of saddle trajectories. To the inventors'knowledge, CB image reconstruction on R-lines of general sourcetrajectories was first suggested by Palamodov, who developed athree-term convolution-based FBP formula, see, e.g., V. P. Palamodov,“Reconstruction from Ray Integrals with Sources on a Curve,” Inv. Prob.,Vol. 20, pp. 239-42, 2004. However, the methods of the present inventionrequire less data than Palamodov's for reconstruction on an R-line andcan accommodate a certain degree of transverse data truncation (inaddition to axial truncation) for reconstruction on some surfaces ofR-lines in the case of the helix or a saddle trajectory.

Another embodiment of a method to compute ƒ on some of the measuredlines that have only one intersection with the source trajectory isdisclosed according to the present invention. This method offers moreflexibility in the design of reconstruction algorithms. In the HCBTcontext, the method enables transverse truncation to be handled over amuch larger region than that allowed by the R-lines approach and alsooffers a way to incorporate redundant data into the reconstructionprocess for possible reduction of image noise.

Additional features and advantages of the invention will be apparentfrom the detailed description which follows, taken in conjunction withthe accompanying drawings, which together illustrate, by way of example,features of embodiments of the present invention.

BRIEF DESCRIPTION OF THE DRAWINGS

The following drawings illustrate exemplary embodiments for carrying outthe invention. Like reference numerals refer to like parts in differentviews or embodiments of the present invention in the drawings.

FIG. 1 is an illustration of the line L(θ,s) and the point r(t,θ,s) thatis on the line L(θ,s) at signed distance t from s according to anembodiment of the present invention.

FIG. 2. is a diagram illustrating the data acquisition with a sourcetrajectory consisting of 3 smooth curves (Γ₁, Γ₂, and Γ₃) according toan embodiment of the present invention.

FIG. 3 is a diagram illustrating projection geometry for a well-orienteddetector according to an embodiment of the present invention.

FIG. 4 is an illustration of the link between the DBP and the Hilberttransform according to an embodiment of the present invention.

FIG. 5 illustrates representations of a detector at three sourcepositions between the endpoints of a central π-line L on a helixaccording to embodiments of the present invention.

FIGS. 6A and 6B illustrate axial views of a single-turn helical scan ofa patient that extends beyond the FOV (dashed circle) in the x-directionaccording to embodiments of the present invention.

FIG. 7A is an exaggerated-pitch illustration of a surface of π-lines anda surface of R-lines (that are not π-lines) on which helicalreconstructions were achieved according to an embodiment of the presentinvention.

FIG. 7B illustrates true images of the FORBILD head phantom without earson the π-line (left side) and the R-line (right side) surfaces accordingto an embodiment of the present invention.

FIG. 8. illustrates reconstructions of the FORBILD head phantom withoutears on the π-line surface of FIG. 7A, according to an embodiment of thepresent invention.

FIG. 9. illustrates reconstructions of the FORBILD head phantom withoutears on the surface of R-lines (that are not π-lines) of FIG. 7Aaccording to embodiments of the present invention.

FIG. 10 illustrates three surfaces of R-lines for the saddle trajectoryaccording to an embodiment of the present invention.

FIG. 11 illustrates reconstructions of the FORBILD head phantom withoutears using embodiments of a method according to the present invention.

FIG. 12 illustrates a wedge volume where accurate reconstruction ispossible over the FOV despite transverse truncation according to anembodiment of the present invention.

FIG. 13 is images representing the Popeye phantom on the central M-linesurface within a wedge according to an embodiment of the presentinvention.

FIG. 14 is images representing the Popeye phantom on a central verticalslice within a wedge as illustrated in FIG. 12 according to anembodiment of the present invention.

FIG. 15 illustrates various reconstructions of the 3-D Popeye phantomdemonstrating the utilization of minimally redundant data for noisereduction according to an embodiment of the present invention.

FIGS. 16A-D are illustrations of the parameters involved in the extendedDBP formula according to embodiments of the present invention.

FIG. 17 illustrates reconstructions from simulated data using theextended DBP equation (52) without (left) and with (right) noise on asurface of R-lines according to embodiments of the present invention.

FIG. 18 is a flow chart of an embodiment of a method of CBreconstruction using truncated CB projections according to the presentinvention.

FIG. 19 is a block diagram of a computer readable media includingcomputer readable instructions configured to implement methods accordingto embodiments of the present invention.

FIG. 20 is a flow chart of a method of image three-dimensional (3D)reconstruction from truncated cone-beam (CB) projections according to anembodiment of the present invention.

FIG. 21 is a block diagram of computed tomography (CT) scanner accordingto an embodiment of the present invention.

DETAILED DESCRIPTION

This detailed description describes a flexible new methodology foraccurate cone-beam reconstruction with source positions on a curve (orset of curves). Embodiments of the methods disclosed herein includeinversion formulas that are based on first backprojecting a simplederivative in the projection space and then applying a Hilbert transforminversion in the image space. The local nature of the projection spacefiltering distinguishes this approach from conventionalfiltered-backprojection methods. This characteristic together with adegree of flexibility in choosing the direction of the Hilbert transformused for inversion offers two important features for the design of dataacquisition geometries, reconstruction algorithms and apparatusesemploying the methods of the present invention. First, the size of thedetector necessary to acquire sufficient data for accuratereconstruction of a given region is often smaller than that required bypreviously documented approaches. In other words, more data truncationis allowed. Second, redundant data can be incorporated for the purposeof noise reduction. The validity of the inversion formulas along withthe application of these two properties are illustrated withreconstructions from computer simulated data as disclosed herein. Inparticular, with reference to the helical cone beam geometry context, itis shown that 1) intermittent transaxial truncation has no effect on thereconstruction in a central region, which means that wider patients canbe accommodated on existing scanners and, more importantly, thatradiation expose can be reduced for region of interest imaging and, 2)at maximum pitch the data outside the Tam-Danielsson window can be usedto reduce image noise and thereby improve dose utilization. Furthermore,the degree of axial truncation tolerated by our approach for saddletrajectories is shown to be larger than that of previous methods.

The following disclosure defines the Hilbert transform of a 3-D objectdensity function on a line in space, and discusses how this transformcan be inverted while being only known on the segment of the line whereƒ is nonzero. As will be seen in the next sections, the definitions andresults presented here provide the foundations for a versatile theoryfor image reconstruction from CB projections.

FIG. 1 is an illustration of the line L(θ,s) and the point r(t,θ,s) thatis on it at signed distance t from s. Throughout this disclosure, thedensity function to be reconstructed is either denoted as ƒ, ƒ(x,y,z),or ƒ(x) with x=(x,y,z). Also, the support of ƒ, which is the smallestclosed set enclosing the region where ƒ is nonzero, is denoted Ω. Forall developments, it is assumed that Ω is compact (closed and bounded).Furthermore, it is assumed that ƒ≧0 like any density function and ƒ iscontinuously differentiable everywhere, which implies that ƒ(x) tends tozero when x tends toward the boundary of Ω.

Let θ be a unit vector and s be a point in space. Then let L(θ,s) be theline of direction θ through s, and letr(t,θ,s)= s+tθ   (1)be a parameterization of the points on this line, with t∈(−∞,+∞), seeFIG. 1. Depending on its position and the location and shape of Ω,L(θ,s) may or may not hit Ω. Also, the intersection of L(θ,s) with Ω maybe more than a single line segment when L(θ,s) hits Ω because θ need notbe convex. However, the convex hull of the intersection of L(θ,s) and Ωis always a single line segment when L(θ,s) hits Ω. This line segmentcorresponds to an interval of t values that is often used in herein andis denoted as region [t_(min)(θ,s),t_(max)(θ,s)], or more simply asregion [t_(min),t_(max)], keeping in mind that t_(min) and t_(max) bothdepend on θ and s. By definition, ƒ is zero at any position r(t,θ,s)defined with t∉(t_(min),t_(max)), and t_(min) and t_(max) are bothfinite numbers since Ω is bounded. Moreover, note that the convex hullof the intersection of L(θ,s) with Ω is in general smaller than theintersection of L(θ,s) with the convex hull of Ω.

Now letk(t,θ,s)=ƒ r(t,θ,s)  (2)be the function that assigns to any given t∈(−∞,+∞) the value of ƒ atlocation r(t,θ,s) on L(θ,s). We define the Hilbert transform of ƒ on theline L(θ,s) as the Hilbert transform of k(t,θ,s) in t, namely$\begin{matrix}{{K( {t,\underset{\_}{\theta},\underset{\_}{s}} )} = {\int_{- \infty}^{\infty}{\frac{1}{\pi( {t - t^{\prime}} )}{k( {t^{\prime},\underset{\_}{\theta}\quad,s} )}{\mathbb{d}t^{\prime}}}}} & (3)\end{matrix}$where the singularity at t′=t is handled in the Cauchy Principal Valuesense.

By definition, the behavior of K(t,θ,s) depends on the position ofL(θ,s) in space. If L(θ,s) has no intersection with Ω, K(t,θ,s)=0everywhere. On the other hand, if L(θ,s) intersects Ω, K(t,θ,s) isnonzero at almost any t, even though k(t,θ,s) is zero-outside regiont∈(t_(min),t_(max)).

Note that the geometrical meaning of t is retained through the Hilberttransform since the operation of equation (3) is shift-invariant. Fromthis observation, the Hilbert transform of ƒ at location x on the lineL(θ,s) has a clear meaning, namely that of K(t,θ,s)|t=(x−s)·θ.Furthermore, since x is on L(θ,s), this expression is independent of sand will be denoted as $\begin{matrix}\begin{matrix}{{K*( {\theta,\underset{\_}{x}} )} = {{{K( {t,\underset{\_}{\theta},\underset{\_}{s}} )}❘t} = {( {\underset{\_}{x} - \underset{\_}{s}} ) \cdot \theta}}} \\{= {{{K( {t,\underset{\_}{\theta},\underset{\_}{s}} )}❘t} = 0.}}\end{matrix} & (4)\end{matrix}$From equations (1)-(4), $\begin{matrix}{{K*( {\underset{\_}{\theta},\underset{\_}{s}} )} = {- {\int_{- \infty}^{\infty}{\frac{1}{\pi\quad t}{f( {\underset{\_}{x} + {t\quad\theta}} )}{\mathbb{d}t}}}}} & (5)\end{matrix}$Also, note from this equation (5) that K*(θ,x) is odd in θ, i.e.,K*(−θ,x)=−K*(θ,x).

The important parts of the above derivation are relations converting CBprojections into values of K*(θ,x). If these values can be computed forpoints x that lie on the same line L(θ,s), a sampling of K(t,θ,s) isobtained asK*(θ,x)| x=s=tθ=K(t,θ,s)  (6)and if that sampling covers the whole region [t_(min),t_(max)],reconstruction of ƒ on L(θ,s) is possible as explained below.

It is well-known that the Hilbert transformation of K(t,θ,s) yields thenegative of k(t,θ,s), i.e., $\begin{matrix}{{k( {t,\underset{\_}{\theta},\underset{\_}{s}} )} = {- {\int_{- \infty}^{\infty}{\frac{1}{\pi( {t - t^{\prime}} )}{K( {t^{\prime},\underset{\_}{\theta},\underset{\_}{s}} )}\quad{\mathbb{d}t^{\prime}}}}}} & (7)\end{matrix}$but this inversion result disregards knowledge on the support ofk(t,θ,s) and consequently requires K(t,θ,s) to be known for all t. Sincek(t,θ,s) is continuous in t and zero outside (t_(min),t_(max)), thereexists a much less demanding inversion formula from the literature onintegral equations. This equation (7) can be written as follows:$\begin{matrix}{{k( {t,\underset{\_}{\theta},\underset{\_}{s}} )} = {{\frac{{C( {\underset{\_}{\theta},\underset{\_}{s}} )} - {k_{1}( {t,\underset{\_}{\theta},\underset{\_}{s}} )}}{w( {t,\underset{\_}{\theta},\underset{\_}{s}} )}{for}\quad t} \in ( {t_{\min},t_{\max}} )}} & (8)\end{matrix}$where C(θ,s) is a constant at fixed (θ,s), $\begin{matrix}{{{w( {t,\underset{\_}{\theta},\underset{\_}{s}} )} = \sqrt{( {t - t_{\min}} )( {t_{\max} - t} )}}{and}} & (9) \\{{k_{1}( {t,\underset{\_}{\theta},\underset{\_}{s}} )} = {\int_{t_{\min}}^{t_{\max}}{\frac{w( {{t^{\prime}\underset{\_}{\theta}},\underset{\_}{s},} )}{\pi( {t - t^{\prime}} )}\quad{K( {t^{\prime},\underset{\_}{\theta},\underset{\_}{s}} )}{{\mathbb{d}t^{\prime}}.}}}} & (10)\end{matrix}$

Equation (8) shows that the portion of K(t,θ,s) defined witht∈(t_(min),t_(max)) allows the recovery of k(t,θ,s) over its wholesupport up to a constant C(θ,s). This recovery essentially involvestaking a simple Hilbert transform of a weighted version of K(t,θ,s) overthe region t∈(t_(min),t_(max)). It is, therefore, as efficient asequation (7), while less demanding on K(t,θ,s). Also, in the limit wheret_(min) and t_(max) tend toward −∞ and +∞, respectively, equation (8)converges to equation (7) for a bounded C(θ,s), since w(t′,θ,s)/w(t,θ,s)converges toward 1. Equation (8) is a generalized version of equation(7).

The presence of the unknown, C(θ,s) in equation (8) shows, however, thatk(t,θ,s) is not uniquely determined by the values of K(t,θ,s) over theregion t∈(t_(min),t_(max)). Additional information is needed to findC(θ,s). One approach is to use the sum of ƒ on the line L(θ,s). That is,if $\begin{matrix}{{g\quad{L( {\underset{\_}{\theta},\underset{\_}{s}} )}} = {\int_{t_{\min}}^{t_{\max}}{{k( {t,\underset{\_}{\theta},\underset{\_}{s}} )}{\mathbb{d}t}}}} & (11)\end{matrix}$is known, then from equation (8) $\begin{matrix}{{C( {\underset{\_}{\theta},\underset{\_}{s}} )} = {\frac{{{gL}( {\underset{\_}{\theta},\underset{\_}{s}} )} + {\int_{t_{\min}}^{t_{\max}}{\frac{k_{1}( {t,\underset{\_}{\theta},\underset{\_}{s}} )}{w( {t,\underset{\_}{\theta},\underset{\_}{s}} )}{\mathbb{d}t}}}}{\int_{t_{\min}}^{t_{\max}}{\frac{1}{w( {t,\underset{\_}{\theta},\underset{\_}{s}} )}{\mathbb{d}t}}}.}} & (12)\end{matrix}$Alternatively, the vanishing of k(t,θ,s) at t=t_(min) and t=t_(max) canbe used to getC(θ,s)=k ₁(t _(min) ,θ,s)=k ₁(t _(max) ,θ,s)  (13)Both approaches are mathematically valid but present some challenges fornumerical implementation. On one hand, equation (12) involves(integrable) singularities at t=t_(min) and t=t_(max) that require somecareful numerical processing. On the other hand, equation (13) requiresa very accurate knowledge of k₁(t,θ,s) at t=t_(min) or t=t_(max), whichis not likely to be available in practice, due to discretization errorsand data noise propagation.

Let ε be a small positive number and let t_(min) ^(ε)=t_(min) ^(ε)−ε andt_(max) ^(ε)=t_(max) ^(ε)+ε. To overcome numerical problems, wesubstitute t_(min) ^(ε) for t_(min) and t_(max) ^(ε) for t_(max) in theright hand side of both equations (9) and (10) to obtain new quantitiescalled w_(ε)(t,θ,s) and k₁ ^(ε)(t,θ,s). This substitution requiresK(t,θ,s) to be known over a slightly larger domain but simplifies thenumerical implementation because equation (8) may be replaced by$\begin{matrix}{{k( {t,\underset{\_}{\theta},\underset{\_}{s}} )} = {{{\chi_{ɛ}(t)}\frac{{C_{ɛ}( {\underset{\_}{\theta},\underset{\_}{s}} )} - {k_{1}^{ɛ}( {t,\underset{\_}{\theta},\underset{\_}{s}} )}}{w_{ɛ}( {t,\underset{\_}{\theta},\underset{\_}{s}} )}{for}\quad t} \in \lbrack {t_{\min}^{ɛ},t_{\max}^{ɛ}} \rbrack}} & (14)\end{matrix}$where C_(ε)(θ,s) is a new constant to be determined and where χ_(ε)(t)is any function that tends to zero when t tends to t_(min) ^(ε) andt_(max) ^(ε) and is equal to one in the region [t_(min),t_(max)] wherek(t,θ,s) may be nonzero. For example $\begin{matrix}{{X_{ɛ}(t)} = \{ \begin{matrix}{1,} & {{{if}\quad{\tau(t)}} < d} \\{\cos^{2}( \frac{\pi( {d - {\tau(t)}} )}{2_{ɛ}} )} & {{{{if}\quad d} < {\tau(t)} < {d + ɛ}},} \\{0,} & {{{if}\quad{\tau(t)}} > {d + ɛ}}\end{matrix} } & (15)\end{matrix}$with d=(t_(min),t_(max))/2 and σ(t)=|t−(t_(min)+t_(max)))/2|. Repeatingthe process that led to equation (12) using equation (14) instead ofequation (8) as a starting point, we get $\begin{matrix}{{C_{\varepsilon}( {\underset{\_}{\theta},\underset{\_}{s}} )} = \frac{{g_{L}( {\underset{\_}{\theta},\underset{\_}{s}} )} + {\int_{t_{\min}^{ɛ}}^{t_{\max}^{ɛ}}{\frac{{\chi_{ɛ}(t)}{k_{1}^{ɛ}( {t,\underset{\_}{\theta},\underset{\_}{s}} )}}{w_{ɛ}( {t,\underset{\_}{\theta},\underset{\_}{s}} )}{\mathbb{d}t}}}}{\int_{t_{\min}^{ɛ}}^{t_{\max}^{ɛ}}{\frac{\chi_{ɛ}(t)}{w_{ɛ}( {t,\underset{\_}{\theta},\underset{\_}{s}} )}{\mathbb{d}t}}}} & (16)\end{matrix}$which does not include any singularity. Or following the ideas that ledto equation (13), we get $\begin{matrix}{{C_{ɛ}( {\underset{\_}{\theta},\underset{\_}{s}} )} = {{\frac{1}{2ɛ}\lbrack {\int_{t_{\min}^{ɛ}}^{t_{\min}}{+ \int_{t_{\max}}^{t_{\max}^{ɛ}}}} \rbrack}{k_{1}^{ɛ}( {t,\underset{\_}{\theta},\underset{\_}{s}} )}{{\mathbb{d}t}.}}} & (17)\end{matrix}$Note that χ_(ε)(t) does not need to be smooth and could in particular bechosen as 1 for τ(t)<d and 0 otherwise, however, for the numericalimplementation it is easier to use a smooth χ_(ε)(t) like that ofequation (15).

The two approaches, equations (16) and (17), have been tested by theinventors and as disclosed in a separate paper on two-dimensional (2-D)classical tomography, see, F. Noo, R. Clackdoyle and J. D. Pack, “ATwo-Step Hilbert Transform Method for 2D Image Reconstruction,” Phys.Med. Biol., Vol. 49, pp. 3903-23, 2004. Both approaches appear toperform equally well when ε is large enough to include a statisticallyrepresentative set of samples of k₁ ^(ε)(t,θ,s) over the regions t_(min)^(ε)<t<t_(min) and t_(max)<t<t_(max) ^(ε). Otherwise, equation (16) ismore robust than equation (17). All the simulations disclosed belowemploy equation (16).

The following discussion defines the CB projections from which thereconstruction of ƒ is to be achieved, then introduces the concept of adifferentiated backprojection (DBP) and shows how this concept links theCB projections to the Hilbert transform of ƒ on measured lines.

Let Ω⁺ be a bounded and convex neighborhood of Ω. The results disclosedherein apply to CB projections measured on a source trajectory thatconsists of a union of N≧1 smooth curves Γ_(l),l=1, . . . ,N such that aneighborhood of each curve lies outside Ω⁺. To simplify the notation,this trajectory is described using a single parameter λ. The sourceposition at λ is a(λ) and the domain of λ is a union of N disjointintervals Λ_(l),l=1, . . . ,N each of which corresponds to one of thecurves Γ_(l) composing the trajectory, see, e.g., FIG. 2.

FIG. 2. is a diagram illustrating the data acquisition with a sourcetrajectory consisting of 3 smooth curves (Γ₁, Γ₂, and Γ₃). The sourcetrajectory is described using a single parameter, Γ_(l) and the sourceposition at Γ_(l) is a(λ). Also, the volume where ƒ is nonzero is Ω, andthere is a bounded and convex neighborhood, Ω⁺, of Ω, such that aneighborhood of each source trajectory curve lies outside Ω⁺. Each CBprojection can be fully described using a flat detector that interceptsall lines passing through the object. The detector plane is orthogonalto the unit vector, e _(w). The detector elements have coordinates u andυ along unit orthogonal vectors, e _(u) and e _(υ) with (u, υ)=(0, 0) atthe orthogonal projection of a(λ) onto the detector plane.

The smoothness of r, is related to the behavior of the tangent vectora′(λ)=da(λ)/d(λ) i.e., it is assumed that a′(λ) is bounded, continuousand nonzero over the interior of Λ_(l) for any l. The nonvanishingcondition on a′(λ) prohibits corners (angular points) on the Γ_(l)curves. However, the source trajectory itself can have corners, forexample a zigzag trajectory is admissible as it can be built byconnections of line segments which are each smooth curves.Geometrically, the curves Γ_(l) may overlap, intersect each other, or besimply connected.

One simple example trajectory is the circle-plus-line trajectory (see,e.g., F. Noo, R. Clackdoyle, and M. Defrise, “Direct Reconstruction ofCone-Beam Data Acquired with a Vertex Path Containing a Circle,” IEEETrans, Image Process., Vol. 7, No. 6, pp. 854-67, Jun. 1998,) whichconsists of a circle and a line orthogonal to the plane of the circle.This trajectory may be described using N=2, Λ_(l)=[0, 2π) and Λ₂=[2π,4π], with $\begin{matrix}{{\underset{\_}{a}(\lambda)} = \{ \begin{matrix}{ {{R_{1}\cos\quad\lambda},{R_{1}\sin\quad\lambda},0} ),} & {{{if}\quad\lambda} \in \Lambda_{1}} \\{R_{2},0,{( {\lambda - {3\pi}} )\frac{H}{\pi}},} & {{{if}\quad\lambda} \in \Lambda_{2}}\end{matrix} } & (18)\end{matrix}$where R₁, R₂, and H are free parameters (except for the condition thatthe source trajectory must be outside Ω⁺).

The CB projection at position λ is the set of integrals $\begin{matrix}{{g( {\lambda,\underset{\_}{\alpha}} )}{\int_{0}^{\infty}{{f( {{\underset{\_}{a}(\lambda)} + {t\quad\underset{\_}{\alpha}}} )}{\mathbb{d}t}}}} & (19)\end{matrix}$obtained for all possible unit vectors α. However, since α(λ) is outsideΩ⁺ and Ω⁺ is convex, there exists a half sphere of unit vectors α thatdo not point toward the object and for which g is consequently zero.Therefore, the CB projection can be fully described using a flat areadetector placed on the opposite side of the object relative to thesource, in a plane that intercepts all lines that diverge from thesource and go through the object.

When using a flat area detector, the CB projection is denotedg_(d)(λ,u,υ), where u and υ are Cartesian coordinates in the detectorplane. For convenience, the orthogonal projection of the source ontothat plane is chose as the origin (u, υ)=(0, 0). This choice gives$\begin{matrix}{{g_{d}( {\lambda,u,\upsilon} )} = {g( {\lambda,{{\underset{\_}{\alpha}( {\lambda,u,\upsilon} )}{with}}} }} & (20) \\{{\underset{\_}{\alpha}( {\lambda,u,\upsilon} )} = {\frac{1}{\sqrt{u^{2} + \upsilon^{2} + D^{2}}}( {{u\quad{\underset{\_}{e}}_{u}} + {\upsilon\quad{\underset{\_}{e}}_{\upsilon}} - {D\quad{\underset{\_}{e}}_{w}}} )}} & (21)\end{matrix}$where e _(w) is the unit vector orthogonal to the detector plane andpointing toward the source, e _(u) and e _(υ) are orthogonal unitvectors in the direction along which u and υ are measured, respectively,and D is the distance from the source to the detector plane (see FIG.2). In general, e _(u), e _(υ), e _(w) and D depend on λ, although thisdependence is not written explicitly.

There exists an infinite number of possible orientations for a flatdetector. When this orientation is such that e _(u) and e _(w) are,respectively, parallel and orthogonal to a′(λ), we say that the detectoris well-oriented. Note that the requirement that the detector interceptsall lines diverging from the source toward the object is incompatiblewith having the detector well-oriented when the line of direction a′(λ)through a(λ) intersects the object. Hereafter, a reference to awell-oriented detector will always assume tacitly that the line ofdirection a′(λ) through a(λ) does not intersect Ω⁺, which is more thatwhat is needed to avoid ambiguities but is convenient to simplify thedevelopments.

By definition, the CB projection is truncated for a given sourceposition a(λ) whenever there are values of u and υ for whichg_(d)(λ,u,υ) is known (measured) on all lines that diverge from a(λ) andpass through a neighborhood of that set.

The inventive method that allows conversion of CB projections intoHilbert transforms of ƒ on measured lines is a version of afiltered-backprojection procedure, referred to herein as thedifferentiated backprojection (DBP). This method may be applied to anysegment of the smooth curves forming the source trajectory. The methodmay be best understood as a 3-D extension of the 2-D DBP introduced inF. Noo, R. Clackdoyle and J. D. Pack, “A Two-Step Hilbert TransformMethod for 2D Image Reconstruction,” Phys. Med. Biol., Vol. 49, pp.3903-23, 2004. For a well-oriented detector, a DBP appears equivalent tothe application of the filtered backprojection steps of Feldkamp et al.,upon removal of the Hilbert kernel from the ramp filter. Thisessentially corresponds to replacing the ramp filter of Feldkamp et al.,by a differentiation filter. More specifically, a DBP for awell-oriented detector is achieved using the following three methodsteps. See FIGS. 2 and 3 for an illustration of some of the involvedquantities.

-   Step 1) Weight each projection to obtain $\begin{matrix}    {{\overset{\_}{g}( {\lambda,u,\upsilon} )} = {\frac{D}{\sqrt{u^{2} + \upsilon^{2} + D^{2}}}{g_{d}( {\lambda,u,\upsilon} )}}} & (22)    \end{matrix}$-   Step 2) Differentiate each weighted projection in u to obtain    $\begin{matrix}    {{g_{F}( {\lambda,u,\upsilon} )} = {\frac{\partial\quad}{\partial u}{{\overset{\_}{g}( {\lambda,u,\upsilon} )}.}}} & (23)    \end{matrix}$-   Step 3) Apply the weighted CB backprojection step of Feldkamp et al.    to g_(F) over any segment [λ₁, λ₂] of one of the smooth curves Γ_(l)    forming the source trajectory to obtain $\begin{matrix}    {{b( {\underset{\_}{x},\lambda_{1},\lambda_{2}} )} = {\int_{\lambda_{1}}^{\lambda_{2}}{\frac{D{{{\underset{\_}{a}}^{\prime}(\lambda)}}{g_{F}( {\lambda,{u^{*}( {\lambda,\underset{\_}{x}} )},{\upsilon^{*}( {\lambda,\underset{\_}{x}} )}} )}}{\lbrack {( {{\underset{\_}{a}(\lambda)} - \underset{\_}{x}} ) \cdot {\underset{\_}{e}}_{w}} \rbrack^{2}}{\mathbb{d}\lambda}}}} & (24)    \end{matrix}$    at any x∈Ω⁺, where u*(λ,x) and υ*(λ,x) are the detector coordinates    of the line that connects x to a(λ), i.e., $\begin{matrix}    {{u^{*}( {\lambda,\underset{\_}{x}} )} = {{- D}\frac{( {\underset{\_}{x} - {\underset{\_}{a}(\lambda)}} ) \cdot {\underset{\_}{e}}_{u}}{( {\underset{\_}{x} - {\underset{\_}{a}(\lambda)}} ) \cdot {\underset{\_}{e}}_{w}}}} & (25) \\    {{\upsilon^{*}( {\lambda,\underset{\_}{x}} )} = {{- D}\frac{( {\underset{\_}{x} - {\underset{\_}{a}(\lambda)}} ) \cdot {\underset{\_}{e}}_{\upsilon}}{( {\underset{\_}{x} - {\underset{\_}{a}(\lambda)}} ) \cdot {\underset{\_}{e}}_{w}}}} & (26)    \end{matrix}$

FIG. 3 is a diagram illustrating projection geometry for a well-orienteddetector according to an embodiment of the present invention. Point A isthe projection of x and has, thus, coordinates u_(A)=u*(λ,x) andυ_(A)=υ*(λ, _) given by equations (25) and (26). When translating x inthe direction of a′(λ) by h: (1) the projection of x moves from A to Bin the direction of e _(u), yielding υ_(B)=υ_(A) for the υ-coordinate ofB and (2) the projected distance between a(λ) and x along the axis e_(w) remains constant.

By definition, the outcome b of a DBP depends on the λ-interval overwhich the backprojection is carried out. This dependence is written downexplicitly in the arguments of b, using the values λ₁ and λ2 that definethe endpoints of the backprojection range. Furthermore, as stated abovein equation (24), there exists a value of l such that [υ₁,λ₂]⊂Λ₁; theDBP is always defined over one smooth segment of the source trajectory.

The DBP concept is a CB generalization of a method that yields theHilbert transform of ƒ from fan-beam projections in 2-D. However, sincethere may exist several ways to have a detector well-oriented, it seemsnatural to ask whether the DBP is really a well-defined(detector-independent) concept. The rest of this section of the detaileddescription proves this assertion and also explains how to extend thedefinition of the DBP to a flat detector that is not well-oriented.

To understand why the DBP of the present invention is independent of theway the detector is well-oriented, a more general expression needs to beadopted for b, namely $\begin{matrix}{{b( {\underset{\_}{x},\lambda_{1},\lambda_{2}} )} = {\int_{\lambda_{1}}^{\lambda_{2}}{\lbrack {\frac{\mathbb{d}\quad}{\mathbb{d}h}{g_{e}( {\lambda,{\underset{\_}{x} + {h\quad{{\underset{\_}{a}}^{\prime}(\lambda)}} - {\underset{\_}{a}(\lambda)}}} )}} \rbrack_{h = 0}{\mathbb{d}\lambda}}}} & (27)\end{matrix}$where g_(ε) is the homogeneous extension of degree-1 of the CBprojection g in α, i.e., $\begin{matrix}{\begin{matrix}{{g_{e}( {\lambda,\underset{\_}{z}} )} = {\frac{1}{\underset{\_}{z}}{g( {\lambda,\frac{\underset{\_}{z}}{\underset{\_}{z}}} )}}} & {{for}\quad{any}} & \underset{\_}{z}\end{matrix} \neq 0} & (28)\end{matrix}$or, from equations (19) and (28) $\begin{matrix}\begin{matrix}{{g_{e}( {\lambda,\underset{\_}{z}} )} = {\int_{0}^{\infty}{{f( {{\underset{\_}{a}\quad(\lambda)} + {t\underset{\_}{z}}} )}{\mathbb{d}t}}}} & {{for}\quad{any}} & {\underset{\_}{z} \neq 0.}\end{matrix} & (29)\end{matrix}$Equation (27) is clearly detector-independent and, like equation (24),defines the DBP as a description of how a backprojection at a given xvaries when x is slightly moved in the direction of a′(λ) during thebackprojection.

To prove that equation (27) is equivalent to equation (24) for awell-oriented detector, note first from equation (28) that$\begin{matrix}\begin{matrix}{{g_{e}( {\lambda,{\underset{\_}{x} - {\underset{\_}{a}(\lambda)}}} )} = \frac{g_{d}( {\lambda,{u^{*}( {\lambda,\underset{\_}{x}} )},{\upsilon^{*}( {\lambda,\underset{\_}{x}} )}} }{{\underset{\_}{x} - {\underset{\_}{a}(\lambda)}}}} \\{= \frac{\overset{\_}{g}( {\lambda,{u^{*}( {\lambda,\underset{\_}{x}} )},{\upsilon^{*}( {\lambda,\underset{\_}{x}} )}} }{{\underset{\_}{( a }(\lambda)} - {\underset{\_}{ x )} \cdot {\underset{\_}{e}}_{w}}}}\end{matrix} & (30)\end{matrix}$where u*(λ,x) and υ*(λ,x) are the quantities of equations (25) and (26),and where {overscore (g)} is the weighted projection of equation (22),see FIG. 3. The substitution of x+ha′(λ) for x in equation (30) givesaccess to a detector-coordinate expression of the integrand in equation(27). In the particular case of a well-oriented detector, both thedenominator in equation (30) and the expression of v*are unaffected bythis substitution because a′(λ) is parallel to e _(u) and orthogonal toe _(w). So, the dependence in h of the value of g_(e) in equation (27)comes only from a difference in values of u*, i.e., $\begin{matrix}{{g_{e}( {\lambda,{\underset{\_}{x} + {h{\underset{\_}{a^{\prime}}(\lambda)}} - {\underset{\_}{a}(\lambda)}}} )} = \frac{\overset{\_}{g}( {\lambda,{u^{*}( {\lambda,{\underset{\_}{x} + {h{\underset{\_}{a^{\prime}}(\lambda)}}}} )},{\upsilon^{*}( {\lambda,\underset{\_}{x}} )}} }{( {{\underset{\_}{a}(\lambda)} - \underset{\_}{x}} ) \cdot {\underset{\_}{e}}_{w}}} & (31)\end{matrix}$see again FIG. 3. The differentiation of this expression with respect toh followed by a setting of h to zero directly yields equation (24) fromequation (27) because from equation (25) $\begin{matrix}{{u^{*}( {\lambda,{\underset{\_}{x} + {h{\underset{\_}{a^{\prime}}(\lambda)}}}} )} = {{u^{*}( {\lambda,\underset{\_}{x}} )} + {h\frac{D{{\underset{\_}{a^{\prime}}(\lambda)}}}{( {{\underset{\_}{a}(\lambda)} - \underset{\_}{x}} ) \cdot {\underset{\_}{e}}_{w}}}}} & (32)\end{matrix}$for any well-oriented detector.

For a detector that is not well-oriented, the DBP is obtained from thesubstitution of x+ha′(λ2) for x in equation (30) that yields adetector-coordinate expression of equation (27). Straightforward butlengthy calculations show that the DBP for an arbitrarily orienteddetector remains that of equation (24) provided the following expressionfor g_(F) is used instead of equation (23) $\begin{matrix}\begin{matrix}{{g_{F}( {\lambda,u,\upsilon} )} = {{{{\underset{\_}{e}}_{T} \cdot ( {{\underset{\_}{e}}_{u} + \frac{u{\underset{\_}{e}}_{w}}{D}} )}\frac{\partial\overset{\_}{g}}{\partial u}} +}} \\{{{{\underset{\_}{e}}_{T} \cdot ( {{\underset{\_}{e}}_{v} + \frac{u{\underset{\_}{e}}_{w}}{D}} )}\frac{\partial\overset{\_}{g}}{\partial\upsilon}} - {{{\underset{\_}{e}}_{T} \cdot {\underset{\_}{e}}_{e}}\frac{\overset{\_}{g}}{D}}}\end{matrix} & (33)\end{matrix}$where e _(T) is the unit vector along a′(λ), i.e., e _(T)=a′(λ)∥a′(λ)∥.FIG. 4. is an illustration of the link (equation (34)) between the DBPand the Hilbert transform. For x∈Ω, the DBP over [λ₁,λ₂] with theaddition of boundary terms yields the difference between two Hilberttransforms of ƒ at x, namely that on the line L₂ in the direction ofω(λ₂,x), and that on the line L₁ in the direction of ω(λ₁,x).

As stated above, the DBP provides a link between CB projections and theHilbert transform of ƒ on measured lines. Using the equations (1)-(6),this link can be written in the following form: for any x∈Ω⁺,1/πb _(a)( x,λ ₁,λ₂)=K*(ω(λ₂ ,x ), x )−K*(ω(λ₁ ,x ), x )  (34)where $\begin{matrix}{{\underset{\_}{\omega}( {\lambda,\underset{\_}{x}} )} = \frac{\underset{\_}{x} - {\underset{\_}{a}(\lambda)}}{{\underset{\_}{x} - {\underset{\_}{a}(\lambda)}}}} & (35)\end{matrix}$and where $\begin{matrix}{{b_{a}( {\underset{\_}{x},\lambda_{1},\lambda_{2}} )} = {{b( {\underset{\_}{x},\lambda_{1},\lambda_{2}} )} + {\sum\limits_{q = 1}^{2}\frac{( {- 1} )^{q}{g( {\lambda_{q},{\underset{\_}{\omega}( {\lambda_{q},\underset{\_}{x}} )}} )}}{{\underset{\_}{x} - {\underset{\_}{a}( \lambda_{q} )}}}}}} & (36)\end{matrix}$is the DBP at x over the interval [λ₁,λ₂] with added boundary terms.

For a geometric interpretation of equation (34), picture the two linesL₁ and L₂ that connect a(λ₁, ) and a(λ₂, ) to x, respectively, as shownin FIG. 4. Then, picture the Hilbert transform of ƒ at x on both lines,L₁ and L₂, in the direction from a(λ₁, ) to x and from a(λ2, ) to x,respectively. These two directions are those given by ω(λ₁,x) andω(λ₂,x). The Formula in equation (34) states that the difference betweenthese two Hilbert transforms at x is proportional to b_(a)(x,λ₁,λ₂),where b_(a)(x,λ₁,λ₂) is obtained by adding two boundary terms to the DBPover the segment of source trajectory joining a(λ₁) to a(λ₂). Note thatthe two boundary terms added to the DBP are just the ratios of the CBmeasurements on L₁ and L₂ to the distances from x to a(λ₁) and a(λ₂),respectively.

In general, the usefulness of equation (34) for CB reconstruction is notat all obvious since equation (34) only gives access to differencesbetween Hilbert-transform values. However, when there exists an R-linethrough x, the situation becomes wholly different, and equation (34) isfound to deliver two fundamental formulas for CB reconstruction. Thesetwo formulas can be stated as follows:

Formula F1: If a(λ₁) and a(λ₂) are two arbitrary positions on one of thesmooth curves forming the source trajectory, then $\begin{matrix}{{K^{*}( {{\underset{\_}{\omega}( {\lambda_{1},\underset{\_}{x},} )}\underset{\_}{x}} )} = {\frac{- 1}{2\pi}{b_{a}( {\underset{\_}{x},\lambda_{1},\lambda_{2}} )}}} & (37)\end{matrix}$for any x∈Ω⁺ that is on the line joining a(λ₁) to a(λ₂).

Formula F2: If a(λ₁), a(λ₂) and a(λ₃) are three arbitrary positions onone of the smooth curves forming the source trajectory, then$\begin{matrix}{{K^{*}( {{\underset{\_}{\omega}( {\lambda_{3},\underset{\_}{x},} )}\underset{\_}{x}} )} = {\frac{- 1}{2\pi}( {{b_{a}( {\underset{\_}{x},\lambda_{3},\lambda_{1}} )} + {b_{a}( {\underset{\_}{x},\lambda_{3},\lambda_{2}} )}} )}} & (38)\end{matrix}$for any x∈Ω⁺ that is on the line joining a(λ₁) to a(λ₂).

Formulas F1 and F2 each give a way to obtain samples of the Hilberttransform of ƒ on measured lines, and provide thereby a basis for CBreconstruction as discussed previously and expanded in the next section.Note, however, that formula F1 is more restrictive than formula F2 as ityields only samples of the Hilbert transform of ƒ on R-lines. In fact,formula F1 is just the particular case of formula F2 defined with λ₃=λ₁,but formula F1 is sufficiently useful to be individually stated.

To obtain formula F1, as previously stated that K* is odd in its firstargument, and note that ω(λ₂,x)=−ω(λ₁,x) for x on the R-line from a(λ₁)to a(λ₂), Hence, under the conditions of formula F1, the right hand sideof equation (34) is −2K*(ω(λ₁,x),x), and equation (34) is equivalent toequation (37).

Formula F2 is obtained by applying equation (34) twice to get1/πb _(a)( x,λ ₃,λ₁)=K*(ω(λ₁ ,x ), x )−K*(ω(λ₃ ,x ), x )  (39)1/πb _(a)( x,λ ₃,λ₂)=K*(ω(λ₂ ,x ), x )−K*(ω(λ₃ ,x ), x )  (40)Since K*(ω(λ₂,x),x)=−K*(ω(λ₁,x),x) under the conditions of formula F2,the addition of these two equations side to side immediately yieldsformula F2.

A proof of equation (34) is now given. Consider the general expressionof equation (27) for the DBP using the notation G_(e)(λ,x) for theintegrand in that expression. For x∈Ω⁺, equation (29) yields$\begin{matrix}\begin{matrix}{{G_{e}( {\lambda,\underset{\_}{x}} )} = \lbrack {\frac{\mathbb{d}}{\mathbb{d}h}{\int_{0}^{\infty}{{f( {{\underset{\_}{a}(\lambda)} + {t( {\underset{\_}{x} + {h{\underset{\_}{a^{\prime}}(\lambda)}} - {\underset{\_}{a}(\lambda)}} )}} )}\quad{\mathbb{d}t}}}} \rbrack_{h = 0}} \\{= \lbrack {\frac{\mathbb{d}}{\mathbb{d}h}{\int_{- \infty}^{\infty}{{f( {{\underset{\_}{a}(\lambda)} + {t( {\underset{\_}{x} + {h{\underset{\_}{a^{\prime}}(\lambda)}} - {\underset{\_}{a}(\lambda)}} )}} )}\quad{\mathbb{d}t}}}} \rbrack_{h = 0}}\end{matrix} & (41)\end{matrix}$because a(λ) is outside Ω⁺ while x+ha′(λ) is in Ω⁺ for h small enoughsince Ω⁺ is open. We insert the derivative with respect to h inside theintegral of equation (41) and use the chain rule to find $\begin{matrix}{{G_{e}( {\lambda,\underset{\_}{x}} )} = {\int_{- \infty}^{\infty}{t{{\underset{\_}{a^{\prime}}(\lambda)} \cdot ( {\overset{->}{\nabla}\quad f} )}( {{\underset{\_}{a}(\lambda)} + {t( {\underset{\_}{x} - {\underset{\_}{a}(\lambda)}} )}} ){{\mathbb{d}t}.}}}} & (42)\end{matrix}$Then, we apply again the chain rule to find $\begin{matrix}{{\frac{\mathbb{d}}{\mathbb{d}\lambda}{f( \underset{\_}{\xi} )}} = {( {1 - t} ){{\underset{\_}{a^{\prime}}(\lambda)} \cdot ( {\overset{->}{\nabla}\quad f} )}( \underset{\_}{\xi} )}} & (42)\end{matrix}$for ξ=a(λ)=t(x−a(λ)) and, thus, from equation (42) $\begin{matrix}\begin{matrix}{{G_{e}( {\lambda,\underset{\_}{x}} )} = {\int_{- \infty}^{\infty}{\frac{t}{1 - t}\frac{\mathbb{d}}{\mathbb{d}\lambda}{f( {{\underset{\_}{a}(\lambda)} + {t( {\underset{\_}{x} - {\underset{\_}{a}(\lambda)}} )}} )}\quad{\mathbb{d}t}}}} \\{= {\frac{\mathbb{d}}{\mathbb{d}\lambda}{\int_{- \infty}^{\infty}{\frac{t}{1 - t}{f( {{\underset{\_}{a}(\lambda)} + {t( {\underset{\_}{x} - {\underset{\_}{a}(\lambda)}} )}} )}\quad{\mathbb{d}t}}}}}\end{matrix} & (44)\end{matrix}$which yields from t/(1−t)=−1+1/(1−t) and equation (29) $\begin{matrix}\begin{matrix}{{G_{\quad e}( {\lambda,\underset{\_}{x}} )} = \frac{\mathbb{d}}{\mathbb{d}\lambda}} \\{\{ {- {g_{e}( {\lambda,{\underset{\_}{x} - {\underset{\_}{a}(\lambda)} + {\int_{- \infty}^{\infty}\frac{t}{1 - t}}}} }} } \\{ {{f( {{\underset{\_}{a}(\lambda)} + {t( {\underset{\_}{x} - {\underset{\_}{a}(\lambda)}} )}} )}\quad{\mathbb{d}t}} \}.}\end{matrix} & (45)\end{matrix}$Next, the following change of variable $\begin{matrix}{{t = {1 + \frac{\tau}{{\underset{\_}{x} - {\underset{\_}{a}(\lambda)}}}}},{\frac{dt}{1 - t} = \frac{d\quad\tau}{\tau}}} & (46)\end{matrix}$is applied along with equations (5), (28), and (35) to get$\begin{matrix}{{G_{e}( {\lambda,\underset{\_}{x}} )} = {\frac{\mathbb{d}}{\mathbb{d}\lambda}{\{ {\frac{- {g( {\lambda,{\underset{\_}{\omega}( {\lambda,\underset{\_}{x}} )}} )}}{{\underset{\_}{x} - {\underset{\_}{a}(\lambda)}}} + {\pi\quad{K^{*}( {{\underset{\_}{\omega}( {\lambda,\underset{\_}{x}} )},\underset{\_}{x}} )}}} \}.}}} & (47)\end{matrix}$Integrating G_(e)(λ,x) in λ according to the expression of the DBP inequation (27), of which G_(e)(λ,x) is the integrand, yields directlyequation (34) upon some rearrangement of the resulting terms. Note thatthe conditions assumed on ƒ, Ω⁺, and the source trajectory justifies themanipulations hereabove, but may not all be necessary; minimumconditions yielding equation (34) were not investigated.

The following discussion investigates the impact of formulas F1 and F2on reconstruction from CB projections on a source trajectory thatconsists of a single smooth curve. The impact of formula F1 is discussedfirst. Next, the added advantages of formula F2 are investigated. Ineach case, results from computer-simulated data are provided in supportof the theory. As in the previous sections, the theory presented hereapplies to general source trajectories. However, the discussion andexamples are mostly focused on a helical trajectory because a goodknowledge about the properties of the source trajectory is required fora proficient application of our results. Such properties have only beendocumented for the helix and the less familiar saddle trajectories, seee.g., J. D. Pack, F. Noo and H. Kudo, “Investigation of SaddleTrajectories for Cardiac CT Imaging in Cone-Beam Geometry,” Phys. Med.Biol., Vol. 49, No. 11, pp. 2317-36.

Note that we are only interested in recovering ƒ at points inside Ω. Forthis reason, Ω⁺ does not appear explicitly in the discussion andequations below. However, the fact that the DBP is defined over Ω⁺ isheavily used because the convex hull of the intersection of a line withΩ does not generally lie in Ω but in Ω⁺ and the DBP over that convexhull is needed for our developments.

The following is a discussion on reconstruction on R-lines. Recall fromabove that formula F1 provides a means to obtain samples of the Hilberttransform of ƒ on R-lines, and recall from Section II-B that ƒ can berecovered on any line L in space when the Hilbert transform of ƒ on L isknown over the convex hull of the intersection of L with Ω. Combiningthese two results together, a procedure is obtained for reconstructionof ƒ on any R-line that intersects Ω. Let a(λ₁) and a(λ₂) be theendpoints of the R-line, and let θ ₁₂ be the unit vector in thedirection from a(λ₁) to a(λ₂). Following the notation of developedabove, the expression for ƒ on the R-line and its Hilbert transform aregiven by the functions k(t,θ ₁₂,a(λ₁)) and K(t,θ ₁₂,a(λ₁)),respectively, see equations (2) and (3). Using this notation thereconstruction procedure according to an embodiment of the presentinvention may be described in the following three steps.

-   1) Determine the region [t_(min),t_(max)] that defines the convex    hull of the intersection of the R-line with Ω.-   2) Use formula F1 to get the Hilbert transform of ƒ on the R-line    for t∈(t_(min) ,t_(max)), i.e., compute $\begin{matrix}    {{K( {t,{\underset{\_}{\theta}}_{12},{\underset{\_}{a}( \lambda_{1} )}} )} = {{- \frac{1}{2\pi}}{b_{a}( {\underset{\_}{r},( {t,{\underset{\_}{\theta}}_{12},{\underset{\_}{a}( \lambda_{1} )}} ),\lambda_{1},\lambda_{2}} }}} & (48)    \end{matrix}$-    from the CB projections for any t∈(t_(min),t_(max)) where r,(t,θ    ₁₂,a(λ₁)) is the position vector on the R-line given by equation    (1).-   3) Apply equation (8) to obtain k(t,θ ₁₂,a(λ₁)) from K(t,θ ₁₂,a(λ₁))    for any t∈(t_(min),t_(max)).    Note that equation (48) is simply a rewriting of equation (37) using    r,(t,θ ₁₂,a(λ₁)) for x following equation (6).

The above reconstruction procedure is similar to the method suggested byZou and Pan for reconstruction on π-lines in HCBT (helical CBtomography) but is not limited to π-lines nor to a helical trajectory.Even so, it only applies to points on R-lines and is, therefore, toorestrictive for some source trajectories since the set of points whereTuy's condition is satisfied is in general larger than the set of pointsthat belong to R-lines. However, this drawback is offset by a lowrequest on detector coverage. To achieve reconstruction on an R-line asdescribed above, we just need to guarantee the DBP (with boundary terms)in equation (48) can be computed at all required t. Using the visibilityconcept described above, this requirement yields the following datacompleteness condition:

-   -   Accurate reconstruction on an R-line is possible whenever the        convex hull of the intersection of the R-line with Ω is visible        while the source travels between the endpoints of the R-line.

This condition requires little compared to what all previously publishedCB reconstruction methods would require for reconstruction at points onan R-line. A good illustration of this feature is obtained in theparticular case of HCBT, where the CB projection of a π-line from asource position between its endpoints is just a line segment connectinga point on the lower boundary of the TD (Tam-Danielsson) detector windowto a point on the upper boundary of this window. See FIG. 5 anddiscussion below which compares at various source positions between theendpoints of a central π-line the detector region required by the aboveprocedure for reconstruction on that π-line. The difference is clearlysignificant. Looking at this difference from a global viewpoint, Zou andPan stated that accurate reconstruction of the entire FOV with the same(low) overscan as Katsevich's formula is possible using only the datainside the TD window. In addition to this global viewpoint, we note fromFIG. 5 that there are projections in which lateral parts of the FOV(and, thus, the imaged object) need not be irradiated for reconstructionon a π-line, i.e., transverse truncation is admissible at some sourcepositions, unlike with Katsevich's formula.

FIG. 5 illustrates representations of a detector at three sourcepositions between the endpoints of a central π-line L on a helixaccording to embodiments of the present invention. The light gray regionis the projection of a central cylinder that encloses Ω. Reconstructionof ƒ on the intersection of L with this cylinder using the R-line methoddisclosed herein requires data in the black region. This black region ismuch smaller than that needed to apply Katsevich's formula, i.e., thedark gray region.

The inventors have closely examined the R-line reconstruction proceduredescribed herein in terms of transverse truncation in HCBT and have madethe following observations. Consider a helical scan with a pitch smallenough for the detector area to include the TD window, and consider apatient that extends outside the FOV in the x-direction (the directionorthogonal to sagittal slices) as illustrated in FIG. 6A. Given theprojection geometry, the DBP is computable at any point in the FOV.Furthermore, the π-lines parallel to the (y,z)-plane through the FOVhave their intersection with the patient almost always completelyincluded in the FOV, as shown in FIG. 6A. Hence, accurate reconstructionis possible on the shaded surface on π-lines seen in that figure, eventhough the projections are transversely truncated on both sides at somesource positions. Unfortunately, the situation is different for π-linesparallel to the (x, z)-plane through the FOV, as illustrated in FIG. 6B.The intersection of most of these π-lines with the patient extendsoutside the FOV, preventing a complete computation of the DBP on thatintersection, and making it impossible to accurately reconstruct ƒanywhere on the π-line according to the method of the present inventionas disclosed in this section. Combining positive and negativeobservations, we can state that despite the presence of transversetruncation in some projections, the R-line reconstruction proceduredisclosed above allows accurate reconstruction over the volume definedby the π-lines whose intersection with the patient is fully included inthe FOV. Depending on the amount of truncation, this volume may or maynot be significant.

FIGS. 6A and 6B illustrate axial views of a single-turn helical scan ofa patient that extends beyond the FOV (dashed circle) in the x-directionaccording to embodiments of the present invention. In FIG. 6A the datais complete for an R-line reconstruction (the method of Section IV-A) onmost π-lines parallel to the (y, z)-plane through the FOV (see the largeshaded area). In FIG. 6B the data is incomplete for an R-linereconstruction on most π-lines that are parallel to the (x, z)-plane.

Experiments from computer-simulated data of the FORBILD head phantomwithout the ears have been performed to test the accuracy of anembodiment of the R-line reconstruction method of the present inventiondescribed above. See Friedrich-Alexander-University Erlangen-Nuimberg,Institute of Medical Physics, Erlangen, Germany (online), available at:www.imp.uni-erlangen.de/forbild/english/results/index.htmfor adescription of the FORBILD head phantom. These experiments used thehelical trajectory of $\begin{matrix}{{\underset{\_}{a}(\lambda)} = ( {{R\quad\cos\quad\lambda},{R\quad\sin\quad\lambda},\frac{P\quad\lambda}{2\pi}} )} & (49)\end{matrix}$and the saddle trajectory of $\begin{matrix}{{{\underset{\_}{a}(\lambda)} = ( {{R\quad\cos\quad{\phi cos}\quad\lambda},{R\quad\cos\quad{\phi sin}\quad\lambda},{R\quad\sin\quad\phi}} )}{{{with}\quad\phi} = {\arctan( {\frac{P}{R}\cos\quad 2\quad\lambda} )}}} & (50)\end{matrix}$where in each case λ is the polar angle in the (x, z)-plane, and P and Rare shaped parameters. Note that equation (5) was referred to as theX-saddle in J. D. Pack, F. Noo and H. Kudo, “Investigation of SaddleTrajectories for Cardiac CT Imaging in Cone-Beam Geometry,” Phys. Med.Biol., Vol. 49, No. 11, pp. 2317-36.

Helical CB projections were simulated with and without Poisson noisecorresponding to an emission of 300,000 photons per ray, using R=57 cmand P=3.24 cm. There were 1160 projections per turn, the detector planewas parallel to the z-axis at distance D=104 cm from the source, theu-axis was parallel to the (x, y)-plane, and the projections werediscretized on a grid of square pixels of size 0.14 cm. Reconstructionswere then achieved on two sets of R-lines parallel to the (x, y)-planeas illustrated in FIG. 7A. The R-lines in the first set were π-lines,while the R-lines in the second set were lines connecting verticesseparated by a polar angle between 2π and 4π. Each set of R-lines formeda surface that was parameterized in (x, y) for display andreconstruction purposes, using square pixels of side 0.075 cm. Thestructures of the phantom over the two surfaces of R-lines are shown inFIG. 7B. Note how the eyes appear to have different sizes as aconsequence of the obliqueness of the R-lines relative to the (x,y)-plane. In an axial (z) slice, the eyes would have the same size.

FIG. 7A is an exaggerated-pitch illustration of a surface of π-lines anda surface of R-lines (that are not π-lines) on which helicalreconstructions were achieved. FIG. 7B illustrates true images of theFORBILD head phantom without ears on the π-line (left side) and theR-line (right side) surfaces. The images are displayed with a compressedgrayscale of 100 HU covering the values from 0 to 100 HU.

FIG. 8. illustrates reconstructions of the FORBILD head phantom withoutears on the π-line surface of FIG. 7A. Reconstructions with (right side)and without (left side) noise added to the data are shown in FIG. 8. Thetop row in FIG. 8 illustrate results obtained using an embodiment of themethod of the present invention. The bottom row in FIG. 8 illustratesresults obtained using the prior art method of Katsevich. The grayscalein FIG. 8 is: 0 to 100 HU.

FIG. 9. illustrates reconstructions of the FORBILD head phantom withoutears on the surface of R-lines (that are not π-lines) of FIG. 7A.Reconstructions with (right side) and without (left side) noise added tothe data are shown. The top row in FIG. 9 illustrates results obtainedusing embodiments of the method of the present invention. The bottom rowin FIG. 9 illustrates results obtained by computing the output ofKatsevich's formula onto the R-line surface. The grayscale in FIG. 9 is:0 to 100 HU.

FIGS. 8 and 9 compare two reconstruction results for the first and thesecond surface of R-lines respectively, obtained using discretizationtechniques similar to those disclosed in F. Noo, J. Pack and D.Heuscher, “Exact Helical Reconstruction Using Native Cone-BeamGeometries,” Plays. Med. Biol., Vol. 48, pp. 3787-818. The result shownin FIG. 8 was computed using the formulas according to methods of thepresent invention. The result shown in FIG. 9 was obtained by computingthe output of Katsevich's formula onto the points forming the surface ofR-lines. In each case, the R-line reconstructions appeared in goodagreement with the theory supporting them. Furthermore, the R-linemethod performed as well as Katsevich's formula on the surface ofπ-lines (FIG. 8), while using less data. For the other surface ofR-lines (FIG. 9), the situation was different, i.e., the R-linereconstruction was noisier while requiring data in the 3π window that isabout 3 times larger than the TD window. For a discussion on the 3πwindow, see R. Proksa, T. Köhler, M. Grass and J. Timmer, “Then-Pi-Method for Helical Cone-Beam CT,” IEEE Trans. Med. Imag., Vol. 19,No. 9, pp. 848-63, September 2000.

While counterintuitive, this increase in noise is not surprising whenone considers what happens as R tends to so while P remains fixed. Inthat limit, the R-lines become identical to the π-lines, and theintervals [π/2, 3π/2], [3π/2,5π/2] and [5π/2,7π/2] of source positionseach provide enough information for computation of the Hilbert transformof ƒ on the R-lines. However, the sign of the Hilbert transform given bythe second interval of data is opposite to that of the other two, sothat the sum of the three Hilbert transforms is a single Hilberttransform. As a consequence, the reconstruction on the R-lines has threetimes more noise power than the reconstruction on the π-lines since thenoise is additive.

FIG. 10 illustrates three surfaces of R-lines for the saddle trajectory.The R-lines forming each of the three surfaces are parallel to a planecontaining the z-axis. The angle between this plane and the x-axis isrespectively 0, 45, and 90 degrees. The z-direction has been exaggeratedin FIG. 10 for display purposes.

Saddle CB projections were simulated with and without Poisson noise. Theadded noise and the simulation parameters were the same as for thehelical data, except that the value of P was 4.4 cm and the normal tothe detector plane was chosen along the line joining the source positionto the origin x=0. From these projections, reconstructions were achievedon 300 surfaces of R-lines, three of which are illustrated in FIG. 10.The R-lines forming any given surface were parallel to a plane passingthrough the z-axis, and the 300 surfaces differed from each other in theangle between that plane and the x-axis, which varied uniformly from 0to π.

FIG. 11 illustrates reconstructions of the FORBILD head phantom withoutears using embodiments of a method according to the present invention.Reconstructions with (right) and without (left) noise added to the dataare shown. The top of FIG. 11 illustrates results on a surface ofR-lines oriented at 45 degrees from the (x, z)-plane. The bottom of FIG.11 illustrates results on the central sagittal slice created byinterpolation of reconstructions on 300 R-line surfaces with angles withthe (x, z)-plane uniformly distributed between 0 and π. The grayscale inFIG. 11 is: 0 to 100 HU. FIG. 11 shows the reconstruction on the surfaceof R-lines at 45 degrees from the x-axis, and shows a reformatted sliceat x=0 obtained on square pixels of side 0.075 cm using interpolationbetween the 300 reconstructions. As in the helical case, the accuracy ofthe results nicely supports the theory disclosed herein. Furthermore, itcan be seen that there is no particular problem arising from thecombination of reconstructions on different surfaces of R-lines. Theaverage z distance between these surfaces was z=0.0577 cm and eachsurface of R-lines was parameterized using rotated (x, y) coordinatesand discretized using square pixels of side 0.075 cm. The expectedoutcome of the reconstructions can be inferred from FIG. 7B. Theinventors calculated that the R-line method produces here a reduction of26% in axial extent of the required detector area, compared to theapplication of the first reconstruction step of the methods suggested inKatsevich or Chen or Pack-Noo-Kudo cited above. All reconstructionsshown in this section assumed Ω to be a circular cylinder of radius12.75 along the z-axis, and were performed using ε=0.2 cm for inversionof the partial Hilbert transforms using the equations disclosed above.

The following is a discussion of reconstruction on measured lines(M-lines). Let E be the subset of points in Ω⁺ that belong to an R-line,and let an M-line be any measured line that intersects Ω and has theconvex hull of its intersection with Ω included in E. By construction,an M-line need not be an R-line, and formula F2 gives a means to obtainsamples of the Hilbert transform of ƒ on any M-line. Combining thismeans with the results disclosed above, a method is obtained toreconstruct ƒ on any M-line. Let a(λ) and α be, respectively, the sourceposition and the unit directional vector defining an M-line, and letk(t,α,a(λ)) and K(t,α,a(λ)) be respectively the expression for ƒ and itsHilbert transform on this line, following equations (2) and (3). Usingthis notation, the reconstruction method is described in the followingfour steps:

-   1) Determine the region [t_(min),t_(max)] that defines the convex    hull of the intersection of the M-line with Ω.-   2) For each t∈(t_(min),t_(max)), find source positions λ*₁ and λ*₂    that define the endpoints of the R-line through point r(t,α,a(λ)) on    the M-line, see equation (1). By definition, λ*₁ and λ*₂ both depend    on t, α and λ, although the dependence is not written explicitly.-   3) Use formula F2 to get the Hilbert transform of ƒ on the M-line    for any t∈(t_(min),t_(max)) i.e., compute $\begin{matrix}    {{K( {t,\underset{\_}{\alpha},{\underset{\_}{a}(\lambda)}} )} = {{- \frac{1}{2\pi}}( {{b_{a}( {{\underset{\_}{r}( {t,\underset{\_}{\alpha},{\underset{\_}{a}(\lambda)}} )},\lambda,\lambda_{1}^{*}} )} + {b_{a}( {{\underset{\_}{r}( {t,\alpha,{\underset{\_}{a}(\lambda)}} )},\lambda,\lambda_{2}^{*}} )}} )}} & (51)    \end{matrix}$-   4) Apply equation (8) to obtain k(t,α,a(λ)) from K(t,α,a(λ)) for any    t∈(t_(min),t_(max)).    Note that equation (51) is simply a rewriting of equation (38) using    λ for λ₃, and using r(t,α,a(λ)) for x following equation (6).

The above method defines a generalization of the R-line approachdiscussed previously. In comparison with the R-line approach, the M-linereconstruction is still limited to points in E, but there is now adegree of freedom in the selection of the line used for reconstructionat a given location x since this line need not be an R-line. Because thelast reconstruction step is nonlocal, this degree of freedom increasesimaging capabilities in the presence of truncation. For reconstructionon an M-line, we have the following data completeness condition:

-   -   Accurate reconstruction on an M-line through a given source        position a is possible whenever each point x on the convex hull        of the intersection of the M-line with Ω is visible while the        source travels from a to both endpoints of an R-line through x.        This condition differs from that for the R-line approach in that        reconstruction at a given x no longer requires reconstruction to        be possible at every point of the convex hull of the        intersection of Ω and the R-line through x. Thus, flexibility is        added.

Consider again the truncation problem illustrated in FIG. 6 for HCBT.The inventors now that the added flexibility of the M-line method allowsaccurate reconstruction to be achieved in most of the FOV, which is asubset of E by property of the π-lines, see Danielsson et al. Let x beany point that has an (x, y) position within the shaded region of FIG.6A. Geometrically, there always exists a π-line through this point andone source position between the endpoints of this π-line that has thesame x-coordinate as x. Let a(λ*) be that source position. Due to thelocation of a(λ*), the M-line that goes through x from a(λ*) has theconvex hull of its intersection with Ω entirely included in the FOV.Furthermore, this M-line intersects the detector within the TD window.Hence, data within the transversely truncated TD window that covers theFOV of FIG. 6A is sufficient to reconstruct ƒ at any point that has an(x, y) position within the shaded region of FIG. 6A.

D. M-Line Simulations

The inventors have tested the above HCBT result on transverse truncationusing computer simulated data of an intermittently truncated helical CBscan of an abdomen phantom described in Table I. TABLE I POPEYE PHANTOMDESCRIPTION Type x_(c) y_(c) z_(c) a b c θ Density (HU) C 0 0 0 16.5 1010 0 0 C −21 0 0 4.5 7 10 0 0 C 21 0 0 4.5 7 10 0 0 C −21 0 0 2 2 10 0500 C −21 0 0 1.4 1.4 10 0 20 C 21 0 0 2 2 10 0 500 C 21 0 0 1.4 1.4 100 20 E −8.5 −0.5 0 6 5 3.5 −20 60 E −2.5 1.5 1 5 4 3.5 0 60 E −3 1.5 04.3 0.9 1.5 0 55 E 9 0.5 0 6 4.5 5 −30 50 E 8 0.5 0 4 3.2 3.5 −30 30 E0.8 8 −2.5 1.8 1.3 0.5 0 40 E 0.8 7 0 1.8 1.3 0.5 0 40 E 0.8 8 2.5 1.81.3 0.5 0 40 C 0 −7.2 −2 2 1.3 0.9 0 500 C 0 −7.2 −2 1 1 0.9 0 20 C 0−7.2 −0.8 1 1 0.3 0 50 C 0 −7.2 0.4 2 1.3 0.9 0 500 C 0 −7.2 0.4 1 1 0.90 20 C 0 −7.2 1.6 1 1 0.3 0 50 C 0 −7.2 2.8 2 1.3 0.9 0 500 C 0 −7.2 2.81 1 0.9 0 20 C 0 −7.2 4 1 1 0.3 0 50 C 0 −7.2 5.2 2 1.3 0.9 0 500 C 0−7.2 5.2 1 1 0.9 0 20 E 0 −4 0 0.32 0.45 1.5 0 30 E 1 −4 0 0.45 0.32 1.50 30

This phantom, dubbed the Popeye phantom, was specifically designed withlow contrast abdomen structures and large arms to demonstrate theability to handle transverse truncation. All parameters used for theexperiment were identical to those used for the helical simulations ofdiscussed above except that the number of photons used to simulatednoisy data and P were increased to 500,000 and 7.4 cm, respectively.Also, the number of detector elements per row was only large enough tocover a cylindrical FOV of radius 13 cm while the radius needed to avoidtransverse truncation with the Popeye phantom is 25.5 cm.

FIG. 12 illustrates a wedge volume where accurate reconstruction ispossible over the FOV despite transverse truncation according to anembodiment of the present invention. Surfaces A and C are surfaces ofπ-lines, while surface B is a surface of M-lines that are not R-lines.The FOV can be covered using a stack of edges such as that displayedhere. Hence, reconstruction from transversely truncated data is possibleat any z in the FOV. Reconstruction was achieved on a wedge composed ofM-line surfaces three of which are illustrated (though not to scale) inFIG. 12. Stacking several such wedges with alternating orientation (+y,−y, +y . . . ) allows full coverage of the FOV in the axial direction,see e.g., Tam et al. and Danielsson et al. The wedge was discretizedinto 55 M-line surfaces such that the axial distance between the centralpoints on adjacent surfaces was 0.0685 cm. Notice that the top andbottom M-line surfaces (labeled A and C in FIG. 12) are π-line surfacesbut the other 53 including the central one labeled B are not.

FIG. 13 is images representing the Popeye phantom on the central M-linesurface within a wedge, see e.g., surface B in FIG. 12. The top rowshown in FIG. 13 includes (left) original, (right) Katsevichreconstruction with transversely truncated data. The middle row shown inFIG. 13 includes (left) reconstruction from transversely truncated datausing the Hilbert transform approach, (right) Katsevich reconstructionfrom nontruncated data. The bottom row shown in FIG. 13 includes same asmiddle row but with Poisson noise added to the data. The grayscale inFIG. 13 is −80 to 140 HU.

FIG. 14 is images representing the Popeye phantom on a central verticalslice within a wedge as illustrated in FIG. 12. Top row: (left)original, (right) Katsevich reconstruction with transversely truncateddata. Middle row: (left) reconstruction from transversely truncated datausing the Hilbert transform approach, (right) Katsevich reconstructionfrom nontruncated data. Bottom row: same as middle row but with Poissonnoise added to the data. The grayscale in FIG. 14 is −40 to 70 HU.

FIGS. 13 and 14 show a reconstruction on this central slice and on acentral vertical slice also visible in FIG. 12, using a compressedgrayscale (−40 to 70 HU). In both FIGS. 13 and 14, the upper left imageis the original phantom (shown for reference), the upper right is theresult of naively applying the prior art formula of Katsevich despitethe presence of transverse truncation in the data, and the other tworows or images compare the results of the M-line method using truncateddata (left) to those of the Katsevich formula using nontruncated data(right), without and with Poisson noise added to the data. Despiteintermittent transverse truncation of the data, the images producedusing the M-line method are in good agreement with those produced byapplying the exact formula of Katsevich to nontruncated data. Thisobservation is striking when one considers the implications for clinicalCT that accompany it. Specifically, if we assume that the region ofinterest (ROI) can be appropriately positioned, radiation can be morenarrowly collimated in the transverse direction leading to a reductionin the X-ray dose and in the likelihood that the patient will contractcancer as a result of the scan.

The following is another important application of the M-line method toHCBT. Given the new freedom that the M-line method affords with regardto the direction along which the Hilbert inversion is performed forreconstruction at a given location, it is natural to wonder what effect,if any, this choice has on the noise properties of the finalreconstruction. In the case of a 180° 2-D parallel-beam scan, nosignificant difference is observed when the inversion is performed alonglines parallel to the x-axis in stead of the y-axis, see e.g.,Noo-Clackdoyle-Pack cited above. However, in HCBT with athird-generation multi-row detector the situation is different becauseredundant data is always available at pitch values that provide completedata, i.e., data outside the TD window is always measured, see e.g.,Proksa et al. cited above and D. Heusher, K. Brown and F. Noo,“Redundant Data and Exact Helical Cone-Beam Reconstruction,” Phys, Med.Biol., Vol. 49, 2004. Such redundant data can be incorporated into thereconstruction by applying the M-line method to M-lines that falloutside the TD window.

Let λ_(in) and λ_(out) be the first and the last source position atwhich a given point x is visible. Assuming no transaxial truncation anda FOV radius and pitch P low enough to avoid interrupted illumination, xis visible for all λ∈[λ_(in),λ_(out)] and assuming the data is complete,the endpoints λ₁ and λ₂ of the π-line through x are such thatλ_(in)<λ₁<λ₂<λ_(out). Hence, reconstruction at x with the M-line methodcan be achieved using any M-line through x with λ∈[λ_(in),λ_(out)], andin particular λ=λ_(in) and λ=λ_(out), without introducing CB artifacts.Such nonapproximate incorporation of an arbitrary amount of redundantdata in the reconstruction represents a novel achievement, but is onlyof value if it can improve noise performance. The fact that the use ofredundant data in the π-window using the R-line method resulted in anincreased noise level (FIG. 9) suggests that reconstructions usingM-lines outside the TD window might have the same negative result.Experimental results shown in FIG. 15 confirm this suspicion. However,as discussed below, the average of reconstructions obtained fromdifferent M-lines can lead to a significant improvement in noiseperformance thanks to statistical differences that originate from thedifference in the data involved for each reconstruction.

FIG. 15 illustrates various reconstructions of the 3-D Popeye phantomdemonstrating the utilization of minimally redundant data for noisereduction. The grayscale shown in FIG. 15 is −80 to 140 HU. The imagesin FIG. 15 represents the Popeye phantom on a (limited) set of π-linesparallel to the (x, z)-plane and similar to those of FIG. 6B. The topimage in FIG. 15 is the original phantom while the remaining six imagesare reconstructions from simulated data with Poisson noise. Allsimulation and reconstruction parameters are identical to those used forthe HCBT experiments discussed above except that the number of photonswas increased to 500,000 and Ω was assumed to be a central cylinder ofradius 26.4 cm. Also, the number of simulated detector rows was suchthat P was the maximum pitch for data completeness, which yields a dataredundancy of about 30%, see Heuscher et al. The second through fourthimages show the results obtained when the reconstruction at each point xis done by respectively applying the M-line method to three linesthrough x, namely the π-line, the M-line through a(λ_(in)) and theM-line through a(λ_(out)). The fifth image is the average of theprevious three images, the sixth is obtained by applying the formula ofKatsevich, and the seventh is the average of images three, four and six.Table II, below, summarizes the total relative error in each image aswell as that in each arm and in the central section. TABLE II RELATIVEERRORS ON RECONSTRUCTIONS OF THE POPEYE PHANTOM Region Full Left-armRight-arm Central π-rcn 6.7 7.8 7.6 6.0 π⁻-rcn 10.4 14.0 8.5 9.6 π⁺-rcn10.5 8.9 14.0 9.7 Aver. 1 4.8 5.7 5.7 4.2 Kat. rcn 5.6 5.4 5.4 5.8 Aver.2 4.7 5.3 5.5 4.2The results of Table II and the images in FIG. 15 clearly demonstratefor the first time that (minimally) redundant data can be successfullyused to reduce noise in HCBT without introduction of CB artifacts due toapproximations in the reconstruction process.

The discussion in the previous section disregards trajectories that areformed from multiple smooth curves. If multiple source curves exist wecan reconstruct using each one separately and merge the results.However, this approach is often too restrictive. For example, using thecircle-plus-line trajectory of equation (18) we would only be able toreconstruct in the plane of the circle, while it is known thatreconstruction over a much larger region is possible. See e.g., Tuy;Smith; and Noo-Clackdoyle-Defrise cited above. The following sectionextends the fundamental DBP equation (34) to a more general form in sucha way that data from different source curves can be combined forreconstruction.

Consider a point x∈Ω⁺ and P intervals of source positions, [λ_(p)⁻,λ_(p) ⁺], p=1, . . . P, such that 1) each interval defines a segmentof one of the N smooth curves forming the source trajectory and 2) foreach p=1, . . . P−1, either a(λ_(p) ⁺)=a(λ_(p+1) ⁻) or the lineconnecting a(λ_(p) ⁺) to a(λ_(p+1) ⁻) passes through x. FIGS. 16A-Dillustrates this situation for three different source trajectories: 1)the circle-plus-line trajectory of equation (18), 2) two disconnectedcoplanar circle arcs, and 3) a “curl” trajectory that consists of ashort-scan circular arc connected to two segments of line orthogonal toit. FIGS. 16A-D are illustrations of the parameters involved in theextended DBP formula according to embodiments of the present invention.More specifically, FIG. 16A illustrates the circle-plus-line trajectorywith x above the circle and P=2. FIG. 16B illustrates the curltrajectory with x above the plane of the short-scan and P=2. FIG. 16Cillustrates two arcs of a single circle with x between them in the sameplane and P=2. FIG. 16D illustrates the curl trajectory with x above theplane of the short-scan and P=3. The extended DBP formula allowsreconstruction at x except for the geometry of FIG. 16C.

Now, let k_(p)=−1 when x is between a(λ_(p) ⁺) and a(λ_(p+1) ⁻), andk_(p)=1 otherwise, for p=1, . . . P−1. Using these numbers and theequation (36) of the DBP with boundary terms, we have the followingextended DBP formula for reconstruction at x. $\begin{matrix}{{\frac{1}{\pi}{\sum\limits_{p = 1}^{P}{q_{p}{b_{a}( {\underset{\_}{x},\lambda_{p}^{-},\lambda_{p}^{+}} )}}}} - {q_{P}{K^{*}( {\underset{\_}{w},( {\lambda_{P}^{+},\underset{\_}{x}} ),\underset{\_}{x}} )}} - {K^{*}( {\underset{\_}{w},( {\lambda_{1}^{-},\underset{\_}{x}} ),\underset{\_}{x}} )}} & (52)\end{matrix}$where q₁=1 and${q_{P} = {{\prod\limits_{r = 1}^{p - 1}{k_{r}\quad{for}\quad p}} > 1}},$and where ω(λ₁ ⁻,x) and ω(λ_(p) ^(°),x) are given by equation (35). Thisformula is easily proved from equation (34). Let LHS and RHS be,respectively, the left-hand and the right hand-sides of equation (52).From equation (34), we get $\begin{matrix}\begin{matrix}{{LHS} = {\sum\limits_{p = 1}^{P}{q_{P}\{ {{K^{*}( {\underset{\_}{w},( {\lambda_{p}^{+},\underset{\_}{x}} ),\underset{\_}{x}} )} - {K^{*}( {\underset{\_}{w},( {\lambda_{p}^{-},\underset{\_}{x}} ),\underset{\_}{x}} )}} \}}}} \\{= {{RHS} + {\sum\limits_{p = 1}^{P - 1}{q_{P}\{ {{K^{*}( {\underset{\_}{w},( {\lambda_{p}^{+},\underset{\_}{x}} ),\underset{\_}{x}} )} - {k_{p}{K^{*}( {\underset{\_}{w},( {\lambda_{p + 1}^{-},\underset{\_}{x}} ),\underset{\_}{x}} \}}}} }}}} \\{= {RHS}}\end{matrix} & (53)\end{matrix}$since K*(w(λ_(p) ⁺),x)=K*(w(λ_(p+1) ⁻ x),x) for k_(p)=1 and K*(w(λ_(p)⁺,x),x)=−K*(w(λp+1 ⁻,x),x) for k_(p)=−1 due to K* being an odd functionin its first argument.

As in the case of the nonextended DBP formula, we find the usefulness ofequation (52) to be not at all obvious, except when x,a(λ₁ ⁻) anda(λ_(p) ⁺) are collinear. In that case, if q_(p)=1 when x is betweena(λ₁ ⁻) and a(λ_(p) ⁺), the right-hand side of equation (52) becomes2⁻q_(p)K*(w(λ_(p) ⁺,x),x), and equation (52) yields, therefore, a way toobtain a sample of the Hilbert transform of ƒ on the line througha(λ_(p) ⁺).

Going back to the examples of FIG. 16, the following illustrativeobservations can be made. In FIG. 16A and 16B, P=2 and q₂=1 while x isbetween a(λ₁ ⁻) and a(λ₂ ⁺), therefore, equation (52) gives access tothe Hilbert transform of F on the line from a(λ₁ ⁻) and a(λ₂ ⁺). In FIG.16C, P=2 but q₂=−1 while x is between a(λ₁ ⁻) and a(λ₂ ⁺), so equation(52) delivers nothing, in agreement with the understanding that thesource trajectory is not complete for any reconstruction at x. In FIG.16D, P=3 and q₃=−1, but x is not between a(λ_(and a(λ) ₃ ⁺), therefore,equation (52) gives access to the Hilbert transform of ƒ on the M-linethat diverges from a(λ₃ ⁺).

Basically, the extended DBP equation (52) allows us to build the Hilberttransform of ƒ on R-lines and M-lines that were not reachable whenconsidering separately each smooth curve forming the source trajectory.In some cases, these lines are reached by combining data fromdisconnected source curves using R-lines to jump from one curve toanother. In other cases, these lines are reached by just noting thatcorners in the segment of source trajectory that connects the endpointsof an R-line are admissible in the DBP operation.

From a data completeness point of view, we note that reconstruction onan R-line or an M-line using equation (52) requires no more than theconvex hull of the intersection of that line with Ω to be visible overthe source trajectory. Hence, for such a reconstruction significantadvantages are obtained in terms of data requirement over otherpreviously published CB reconstruction methods. For example, in FIG.16A, the reconstruction at x does not require data over the planesthrough x that do not intersect the circle, unlike the method inNoo-Clackdoyle-Defrise cited above, while allowing axial truncation. Foranother example, in FIG. 16B, the reconstruction on the R-line from a(λ₁⁻) to a(λ₂ ⁺) is possible from axially truncated projections, which is acompletely original result.

The inventors have applied equation (52) to simulated projections of theFORBILD head phantom (without ears) on the curl trajectory of FIG. 16B.All applicable parameters from the experiments above were used in thisexperiment. The detector was wide enough to prevent transversetruncation, but not tall enough to prevent axial truncation. Theshort-scan arc covered 206 degrees with 696 source positions. Each linesegment covered a distance of 5 cm with 25 source positions.Reconstruction was achieved on a surface of R-lines connecting pointa(λ₁ ⁻) of FIG. 16B (at 5 cm from the circular arc) to a(λ₂ ⁺) with (λ₂⁺) spanning an arc of 52.8 degrees from the endpoint of the circular arcclosest to the particular a(λ₂ ⁺) in FIG. 16B. The distance betweenadjacent samples on each R-line was 0.075 cm and the angle betweenadjacent R-lines was chosen such that they intersected the short scan atpoints with an angular separation of 0.15 radians from the point of viewat the center of the short scan.

FIG. 17 illustrates reconstructions from simulated data using theextended DBP equation (52) without (left) and with (right) noise on asurface of R-lines. Each R-line was like the R-line from (λ₁ ⁻) to (λ₂⁺) in FIG. 16B. The surface was built by keeping (λ₁ ⁻) fixed and moving(λ₂ ⁺) counterclockwise from the extremity of the short-scan. Thegrayscale of FIG. 17 is 0 to 100 HU. The images in FIG. 17 show thereconstruction obtained with and without Poisson noise added to thedata.

FIG. 18 is a flow chart of an embodiment of a method 1800 of CBreconstruction using truncated CB projections according to the presentinvention. Method 1800 may include performing 1802 a differentiatedbackprojection (DBP) on segments of truncated CB projections to obtain aportion of a Hilbert transform of a density function, ƒ, along measuredlines. Method 1800 may further include reconstructing 1804 the densityfunction, ƒ, on the measured lines by inverse Hilbert transformation.According to an embodiment of method 1800, the measured lines may beredundantly measured lines (R-lines). According to another embodiment ofmethod 1800, performing 1802 a DBP on segments of the truncated CBprojections, may include determining a region [t_(min), t_(max)] thatdefines a convex hull of an intersection of an R-line with Ω andobtaining a Hilbert transform of the density function, ƒ, on the R-linefor any t∈(t_(min), t_(max)). According to yet another embodiment ofmethod 1800, obtaining a Hilbert transform of the density function, ƒ,on the R-line for any t∈(t_(min), t_(max)), may include computing:$\begin{matrix}{{K( {t,{\underset{\_}{\theta}}_{12},{\underset{\_}{a}( \lambda_{1} )}} )} = {{- \frac{1}{2\pi}}{b_{a}( {\underset{\_}{r},( {t,{\underset{\_}{\theta}}_{12},{\underset{\_}{a}( \lambda_{1} )}} ),\lambda_{1},\lambda_{2}} }}} & (48)\end{matrix}$from the truncated CB projections for any t∈(t_(min),t_(max)), wherer,(t,θ ₁₂, a(λ)) is the position vector on the R-line given by:r (t,θ,s )= s+tθ.  (1)

According to another embodiment of method 1800, reconstructing 1804 thedensity function, ƒ on the R-lines may include obtaining k(t,θ ₁₂,a(λ₁))from K(t,θ ₁₂, a(λ₁)) for any t∈(t_(min),t_(max)). According to anotherembodiment of method 1800, the measured lines may be M-lines. Accordingto another embodiment of method 1800, performing a DBP on segments ofthe truncated CB projections may include determining a region [t_(min),t_(max)] that defines a convex hull of an intersection of an M-line withΩ. The embodiment of method 1800 may further include obtaining a Hilberttransform of the density function, ƒ, on the M-line for any t∈(t_(min),t_(max)). According to yet another embodiment of method 1800, obtaininga Hilbert transform of the density function, ƒ, on the M-line for anyt∈(t_(min), t_(max)), may include for each t∈(t_(min),t_(max)), findingsource positions λ*₁ and λ*₂ that define endpoints of an R-line throughpoint r(t,α,a(λ)) on the M-line. The method 1800 may further includecomputing: $\begin{matrix}{{K( {t,\underset{\_}{\alpha},{\underset{\_}{a}(\lambda)}} )} = {{- \frac{1}{2\pi}}{( {{b_{a}( {{\underset{\_}{r}( {t,\alpha,{\underset{\_}{a}(\lambda)}} )},\lambda,\lambda_{1}^{*}} )} + {b_{a}( {{\underset{\_}{r}( {t,\alpha,{\underset{\_}{a}(\lambda)}} )},\lambda,\lambda_{2}^{*}} )}} ).}}} & (51)\end{matrix}$According to another embodiment of method 1800, reconstructing thedensity function, ƒ, on the M-lines may include obtaining k(t,α,a(λ))from K(t,α,a(λ)) for any t∈(t_(min),t_(max)).

FIG. 20 is a flow chart of a method 2000 of image 3D reconstruction fromtruncated cone-beam (CB) projections according to an embodiment of thepresent invention. Method 2000 may include measuring 2002 the truncatedCB projections. Method 2000 may further include performing 2004 adifferentiated backprojection (DBP) on segments of the truncated CBprojections to obtain a portion of a Hilbert transform of a densityfunction along measured lines. Method 2000 may further includereconstructing 2006 the density function on the measured lines byinversion of a finite Hilbert transform. Method 2000 may further includeforming 2008 a 3D image from the reconstructed density function.

FIG. 19 is a block diagram of a computer readable media 1900 includingcomputer readable instructions 1902 configured to implement methods1800, 2000 according to embodiments of the present invention.

FIG. 21 is a block diagram of computed tomography (CT) scanner 2100according to an embodiment of the present invention. CT scanner 2100 mayinclude an object support 2102 configured for holding an object 2104being examined at least partially within an examination region 2114. CTscanner 2100 may further include a rotating gantry 2106 surrounding theobject support 2102 and configured for rotation about the examinationregion 2114. CT scanner 2100 may further include a source 2108 ofpenetrating radiation disposed on the rotating gantry 2106 for rotationtherewith, the source 2108 of penetrating radiation emitting acone-shaped beam 2110 of radiation that passes through the examinationregion 2114 as the rotating gantry 2106 rotates. CT scanner 2100 mayfurther include an array 2112 of radiation detectors on the rotatinggantry 2106 configured to receive cone-beam (CB) projections from theradiation emitted from the source 2108 of penetrating radiation after ithas traversed the examination region 2114 the array 2112 furtherconfigured for a preselected data acquisition geometry. CT scanner 2100may further include an image reconstruction processor 2116 forreconstructing images of the object 2104 from CB projections collectedby the array 2112 of radiation detectors, the image reconstructionprocessor 2116 configured to implement method 1800 of CB reconstructionas described herein. CT scanner 2100 may further include a monitor 2118in communication with the image reconstruction processor 2116 forviewing reconstructed images of the object 2104. According to additionalembodiments of CT scanner 2100, the preselected data acquisitiongeometry may be any one of helix, saddle trajectories, ortwo-orthogonal-circle orbit.

In summary new methodology has been disclosed for accurate and stablereconstruction from CB projections obtained on a very general class ofsource trajectories. Compared to other reconstruction approaches, thismethodology exhibits attractive features but also some drawbacks. Itskey advantage is a low demand on detector coverage for reconstruction ofa given ROI, which can be converted into increased imaging capabilitieson a given system or into dose reduction techniques. Another advantageis its provision of an original means to incorporate redundant data intothe reconstruction for the purpose of reducing data noise.Unfortunately, one necessary condition for reconstruction at a point xusing this methodology is that x belong to an R-line. For some sourcetrajectories, there are points that do not satisfy this condition, butdo satisfy Tuy's condition. For example, consider the trajectoryresulting from equation (50) if 2λ is replaced by 3λ in the definitionof φ. The convex hull of this trajectory, which is the region whereTuy's condition is satisfied, includes a central portion of the z-axisof length PR/√{square root over (R²+P²)}. However, only one point on thez-axis belongs to an R-line, the point x=0. Fortunately, there are manydata acquisition geometries for which the region where Tuy's conditionis satisfied is completely covered by R-lines. The list includes thehelix, the newly introduced saddle trajectories, and even thetwo-orthogonal-circle orbit of Tuy.

As disclosed herein, several implementation strategies using our newmethodology are possible for reconstruction with a given dataacquisition geometry. These strategies differ in the choice of R-linesor M-lines over which the inversion is carried out and have to betailored to the data acquisition to obtain an optimal reconstructionalgorithm. Also, these strategies do not in general allow forreconstruction directly onto a 3-D rectilinear grid since the samplesfor the post-backprojection Hilbert inversion have to be on R-lines orM-lines. However, there are widely accepted reconstruction methods (suchas AMPR) that also exhibit this feature, see e.g., T. Flohr, K.Stierstorfer, H. Bruder, J. Simon, A. Polacin and S. Schaller, “ImageReconstruction and Image Quality Evaluation for a 16-Slice CT Scanner,”Med. Phys., Vol 30, No. 5, 2003.

The noise performance of embodiments of the new reconstructionmethodology of the present invention was found to vary widely accordingto the implementation strategy. A comparison of the images in FIG. 9demonstrated that reconstruction on a surface of R-lines that are notπ-lines using an embodiment of the method of the present inventionproduces a higher noise level than that produced using Katsevich'sformula (which is already not optimal in terms of noise). It isimportant to note, however, that these poor noise characteristics arenot intrinsic to the methodology presented by this paper. In fact, thediscussion herein demonstrated that in HCBT, reconstructing each pointusing the Hilbert transform along three different M-lines and averagingthe results can significantly improve noise characteristics byincorporating data outside the TD window. Since current multi-slice CTscanners always have an amount of data in excess of the TD window, thisability constitutes a significant advance for HCBT.

The method of reconstruction on M-lines is remarkable in that it solvesthe intermittent transverse truncation problem in HCBT. This allows asmaller detector area to be used (and consequently a lower radiationdose to be delivered) for ROI imaging. Alternatively, wider patients canbe accommodated on currently available clinical scanners. Note that thismethod of dealing with transverse truncation in HCBT can be as efficientas applying Katsevich's formula.

Finally, the inventors have observed that reconstruction from dataacquired on a saddle trajectory can be achieved with the R-linereconstruction method using a significantly smaller detector than themethods of the prior art. In addition, the inventors have reason tobelieve that the R-line method is also more efficient, particularly inachieving high temporal resolution through the use of data from anappropriate limited segment of the saddle for each voxel.

While the foregoing advantages of the present invention are manifestedin the detailed description and illustrated embodiments of theinvention, a variety of changes can be made to the configuration, designand construction of the invention to achieve those advantages. Hence,reference herein to specific details of the structure and function ofthe present invention is by way of example only and not by way oflimitation.

1. A method of cone-beam (CB) reconstruction using truncated CBprojections, comprising: performing a differentiated backprojection(DBP) on segments of truncated CB projections to obtain a portion of aHilbert transform of a density function, ƒ, along measured lines; andreconstructing the density function, ƒ, on the measured lines byinversion of a finite Hilbert transform.
 2. The method according toclaim 1, wherein the measured lines comprise redundantly measured lines(R-lines).
 3. The method according to claim 2, wherein performing a DBPon segments of the truncated CB projections, comprises: determining aregion [t_(min), t_(max)] that defines a convex hull of an intersectionof an R-line with Ω; and obtaining a Hilbert transform of the densityfunction, ƒ, on the R-line for any t∈(t_(min), t_(max)).
 4. The methodaccording to claim 2, wherein obtaining a Hilbert transform of thedensity function, ƒ, on the R-line for any t∈(t_(min), t_(max)),comprises computing:${K( {t,{\underset{\_}{\theta}}_{12},{\underset{\_}{a}( \lambda_{1} )}} )} = {{- \frac{1}{2\pi}}{b_{a}( {\underset{\_}{r},( {t,{\underset{\_}{\theta}}_{12},{\underset{\_}{a}( \lambda_{1} )}} ),\lambda_{1},\lambda_{2}} }}$from the truncated CB projections for any t∈(t_(min),t_(max)), wherer,(t,θ ₁₂,a(λ₁)) is the position vector on the R-line given by:r (t,θ,s )= s+tθ.
 5. The method according to claim 2, whereinreconstructing the density function, ƒ, on the R-lines comprisesobtaining k(t,θ ₁₂,a(λ₁)) from K(t,θ ₁₂,α(λ₁)) for anyt∈(t_(min),t_(max)).
 6. The method according to claim 1, wherein themeasured lines comprise M-lines.
 7. The method according to claim 6,wherein performing a DBP on segments of the truncated CB projections,comprises: determining a region [t_(min), t_(max)] that defines a convexhull of an intersection of an M-line with Ω; and obtaining a Hilberttransform of the density function, ƒ, on the M-line for any t∈(t_(min),t_(max)).
 8. The method according to claim 6, wherein obtaining aHilbert transform of the density function, ƒ, on the M-line for anyt∈(t_(min), t_(max)), comprises: for each t∈(t_(min),t_(max)), findingsource positions λ*₁ and λ*₂ that define endpoints of an R-line throughpoint r(t,α,a(λ)) on the M-line; and computing:${K( {t,\underset{\_}{\alpha},{\underset{\_}{a}(\lambda)}} )} = {{- \frac{1}{2\pi}}{( {{b_{a}( {{\underset{\_}{r}( {t,\alpha,{\underset{\_}{a}(\lambda)}} )},\lambda,\lambda_{1}^{*}} )} + {b_{a}( {{\underset{\_}{r}( {t,\alpha,{\underset{\_}{a}(\lambda)}} )},\lambda,\lambda_{2}^{*}} )}} ).}}$9. The method according to claim 6, wherein reconstructing the densityfunction, ƒ, on the M-lines comprises obtaining k(t,α,a(λ)) fromK(t,α,a(λ)) for any t∈(t_(min),t_(max)).
 10. Computer readable media,comprising computer readable instructions configured to implement themethod of claim
 1. 11. A method of image three-dimensional (3D)reconstruction from truncated cone-beam (CB) projections, comprising:measuring the truncated CB projections; performing a differentiatedbackprojection (DBP) on segments of the truncated CB projections toobtain a portion of a Hilbert transform of a density function alongmeasured lines; reconstructing the density function on the measuredlines by inverse Hilbert transformation; and forming a 3D image from thereconstructed density function.
 12. Computer readable media, comprisingcomputer readable instructions configured to implement the method ofclaim
 11. 13. A computed tomography (CT) scanner, comprising: an objectsupport configured for holding an object being examined at leastpartially within an examination region; a rotating gantry surroundingthe object support and configured for rotation about the examinationregion; a source of penetrating radiation disposed on the rotatinggantry for rotation therewith, the source of penetrating radiationemitting a cone-shaped beam of radiation that passes through theexamination region as the rotating gantry rotates; an array of radiationdetectors on the rotating gantry configured to receive cone-beam (CB)projections from the radiation emitted from the source of penetratingradiation after it has traversed the examination region the arrayfurther configured for a preselected data acquisition geometry; an imagereconstruction processor for reconstructing images of the object from CBprojections collected by the array of radiation detectors, the imagereconstruction processor configured to implement a method of CBreconstruction, the method comprising: performing a differentiatedbackprojection (DBP) on segments of the truncated CB projections toobtain a portion of a Hilbert transform of a density function alongmeasured lines; reconstructing the density function on the measuredlines by inverse Hilbert transformation; and a monitor in communicationwith the image reconstruction processor for viewing reconstructed imagesof the subject.
 14. The CT scanner according to claim 13, wherein thepreselected data acquisition geometry comprises helix.
 15. The CTscanner according to claim 13, wherein the preselected data acquisitiongeometry comprises saddle trajectories.
 16. The CT scanner according toclaim 13, wherein the preselected data acquisition geometry comprisestwo-orthogonal-circle orbit.